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Abstract 

Previous  research  at  AFIT  has  resulted  in  the  development  of  a  DGPS-aided 
INS-based  precision  landing  system  (PLS)  capable  of  meeting  the  FA  A  precision 
requirements  for  instrument  landings.  The  susceptibility  of  DGPS  transmissions  to 
interference /jamming  and  spoofing  must  be  addressed  before  DGPS  may  be  safely 
used  as  a  major  component  of  such  a  safety-of-flight  critical  navigational  device. 
This  thesis  applies  multiple  model  adaptive  estimation  (MMAE)  techniques  to  the 
problem  of  detecting  and  identifying  interference/jamming  and  spoofing  failures  in 
the  DGPS  signal.  Such  an  MMAE  is  composed  of  a  bank  of  parallel  filters,  each 
hypothesizing  a  different  failure  status,  along  with  an  evaluation  of  the  current  prob¬ 
ability  of  each  hypothesis  being  correct,  to  form  a  probability-weighted  average  out¬ 
put.  Performance  for  a  representative  selection  of  navigation  component  cases  is 
examined. 

For  interference/jamming  failures  represented  as  increased  measurement  noise 
variance,  results  show  that,  because  of  the  good  FDI  performance  using  MMAE,  the 
blended  navigation  performance  is  essentially  that  of  a  single  extended  Kalman  filter 
artificially  informed  of  the  actual  interference  noise  variance.  Standard  MMAE  is 
completely  unable  to  detect  spoofing  failures  (modelled  as  a  bias  or  ramp  offset  signal 
directly  added  to  the  measurement).  This  thesis  shows  the  development  of  a  moving- 
bank  pseudo-residual  MMAE  (PRMMAE)  to  detect  and  identify  spoofing  failures. 
Using  the  PRMMAE  algorithm,  spoofing  failures  are  very  effectively  detected  and 
isolated;  the  resulting  navigation  performance  is  equivalent  to  that  of  an  extended 
Kalman  filter  operating  in  a  no-fail  environment. 
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MMAE  DETECTION  OF  INTERFERENCE/ JAMMING  AND 
SPOOFING  IN  A  DGPS- AIDED  INERTIAL  SYSTEM 


1.  Introduction 

1.1  Background 

Aircraft  are  often  required  by  necessity  or  bad  luck  to  fly  under  adverse  weather 
conditions.  The  purposes  of  military  air  operations  are  often  best  accomplished 
under  such  low-visibility  conditions.  Precise,  safe  landings  are  required  even  when 
lack  of  visibility  would  make  landing  by  a  human  pilot  impossible.  The  availability 
and  integrity  of  a  precision  landing  system  is  critical  for  safety  of  flight  during  low 
visibility  conditions. 

The  Instrument  Landing  System  (ILS)  currently  in  use  works  on  a  relatively 
simple  localizer /glide  slope  method.  Guidance  to  the  runway  is  provided  by  highly 
directional  radio  signals.  The  localizer  provides  the  information  needed  for  lateral 
guidance  on  the  approach  path.  The  glide  slope  signal  provides  the  needed  vertical 
guidance.  Figure  1.1  shows  a  representation  of  the  localizer  and  glide  slope  radiation 
patterns.  The  carrier  frequencies  used  in  the  ILS  are  about  llOMHz  and  330MHz 
for  the  localizer  and  glide  slope  signals,  respectively.  The  frequencies  indicated  on 
the  diagram  describe  the  modulation  frequencies  used.  The  effect  of  the  ILS  signals 
just  described  is  to  create  a  virtual  funnel  which  will  guide  an  approaching  airplane 
to  the  runway,  as  is  depicted  in  Figure  1.2. 

The  FAA  has  established  clear  precision  requirements  for  instrument  landings. 
These  requirements  constitute  the  performance  constraints  for  the  current  ILS  as 
well  as  for  the  replacement  GPS-based  precision  landing  system  (PLS)  examined  in 
this  thesis.  The  FAA  requirements  are  shown  in  Category  III  required  navigational 
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Localizer  (as  seen  from  above) 


Figure  1.1  Localizer/Glideslope  Method 


precision  (at  the  200  foot  “decision  height”)  is  roughly  depicted  by  the  ellipse  (two- 
sigma  values)  shown  in  Table  1.1.  Figure  1.3.  The  accuracies  listed  under  Category 
III  represent  the  required  precision  for  a  zero  visibility  landing,  wherein  the  pilot  does 
not  make  visual  contact  with  the  runway  until  after  touchdown.  Conceivably,  this 
precision  would  also  make  completely  autonomous  landings  possible  for  unmanned 
vehicles  or  if  the  pilot  became  incapacitated.  The  localizer /glide  slope  ILS  currently 
in  use  is  capable  of  providing  Category  III  accuracy.  See  Table  1.1. 

The  current  ILS  has  two  major  drawbacks  which  have  motivated  the  desire 
for  an  operationally  certified  CPS-based  precision  landing  system  (PLS).  First,  the 
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Vertical  Accuracy 

Horizontal  Accuracy 

Category  I 

db4.1m 

±17.1m 

Category  II 

±1.7m 

±5. 2m 

Category  III 

±0.6m 

±4.1m 

Table  1.1  Precision  Landing  System  Two-Sigma  Accuracies  [3] 


Figure  1.3  Category  III  Precision  (Two  Standard  Deviations) 

current  ILS  has  a  very  limited  area  of  usefulness  (see  Figure  1.2).  The  area  of 
guidance  is  restricted  to  a  straight  path  extended  for  several  miles  from  the  end  of 
the  runway.  It  would  greatly  contribute  to  the  safety  and  efficiency  of  air  operations 
if  controllers  were  able  to  know  the  precise  locations  and  courses  of  all  airplanes  in 
the  controlled  area.  The  current  ILS  has  no  capacity  to  contribute,  in  this  manner, 
to  the  safety,  guidance,  or  regulation  of  typical  traffic  patterns  of  aircraft  in  various 
stages  of  landing  or  taking  off  around  an  airfield.  The  second  major  drawback  of  the 
current  ILS  is  the  maintenance  cost  of  the  aging  system.  The  Microwave  Landing 
System  (MLS)  was  originally  thought  to  be  the  replacement  for  the  current  ILS. 
The  MLS  however,  works  on  the  same  directional  radiation  principles  and  so  has  the 
same  limited  coverage  drawback  of  the  current  ILS.  Recently  the  MLS  replacement 
plan  has  been  permanently  discarded  due  to  the  high  cost  of  the  upgrade  and  the 
good  potential  of  a  GPS-based  system  for  precision  landings  [33,37,46]. 

The  global  positioning  system  (GPS)  has  been  demonstrated  to  provide  very 
precise  position  measurements,  especially  when  combined  with  other  sensors  [18,36]. 
Figure  1.4  [34]  shows  the  two-sigma  accuracy  capability  of  various  GPS  implemen- 
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tations  which  use  the  standard  positioning  service  (SPS)  LI  GPS  transmission.  The 
GPS  accuracy  ellipse  shown  does  not  reflect  the  scheduled  removal  of  selective  avail¬ 
ability  (SA).  SA  is  a  random  dither  imposed  on  the  GPS  signals  in  order  to  degrade 
the  SPS  precision  available  to  non-authorized  (non-U. S.  military)  users.  Differential 
GPS  (DGPS)  provides  increased  navigation  precision  using  two  closely  located  GPS 
receivers.  One  of  the  receivers  must  be  at  a  known  location  so  that  the  GPS  solution 
errors  may  be  isolated.  The  error  corrections  are  then  transmitted  to  the  second 
receiver,  which  corrects  its  own  navigation  solution.  It  can  be  seen  that,  although 
neither  basic  GPS  nor  DGPS  gives  the  degree  of  accuracy  required  for  precision 
landings,  the  accuracy  of  carrier  phase  GPS  (CPGPS)  is  more  than  sufficient  for 
aircraft  landings.  At  the  time  of  this  writing,  receiver  technology  has  not  yet  made 
CPGPS  available  during  all  phases  of  flight,  especially  during  highly  dynamic  flight. 
However,  for  the  benign  flight  patterns  typical  in  the  area  of  airfields,  the  accuracy 
of  CPGPS  could  be  fully  utilized  with  current  technology. 

The  GPS  has  several  inherent  strengths  and  weaknesses  which  directly  influ¬ 
ence  its  application  to  the  PLS  problem.  The  continual,  global  availability  (inter¬ 
mittent  at  very  high  latitudes)  of  the  GPS  signal,  along  with  the  accuracy  available 
from  CPGPS,  make  a  GPS-based  ILS  appear  to  be  an  improvement  (in  positioning) 
from  the  current  ILS.  A  GPS-based  PLS  would  overcome  all  of  the  major  limitations 
of  the  current  ILS,  while  exceeding  its  current  performance  specifications.  It  can  be 
seen  in  Figure  1.4  that  GPS  alone  (the  outer  ellipse)  and  even  differential  GPS  do  not 
meet  the  accuracy  requirements  for  a  PLS  (the  central  ellipse),  especially  in  the  ver¬ 
tical  direction.  This  vertical  weakness  of  GPS  leads  to  the  desire  to  include  a  radar 
altimeter  in  the  integrated  PLS  system.  A  ground-based  pseudolite  will  also  help 
with  the  vertical  precision  problem.  A  pseudolite  improves  the  navigation  solution 
by  providing  GPS  satellite-like  signals  from  a  known  position  (surveyed,  eliminating 
satellite  position  uncertainty)  below  the  user’s  horizon. 
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Figure  1.4  GPS  Accuracies  (Two  Standard  Deviations) 

A  significant  weakness  of  GPS,  with  respect  to  the  PLS  application,  is  the 
susceptibility  of  the  very  low  power  GPS  signal  to  interference.  For  military  applica¬ 
tions  this  interference  could  originate  from  enemy  interference/ jamming  or  spoofing. 
Benign  interference  such  as  that  from  unregulated  electronic  devices  would  proba¬ 
bly  have  the  greatest  effect  on  civil  aviation,  although  interference  from  low-wattage 
jammers  placed  by  terrorists  is  a  viable  concern.  The  susceptibility  of  the  GPS  signal 
to  external  interference  strongly  motivates  the  use  of  a  GPS-aided  INS-based  PLS, 
rather  than  a  standalone  GPS-based  PLS.  The  navigation  solution  of  a  low-quality 
INS  drifts  at  a  rate  of  about  four  nautical  miles  per  hour  or  400  feet  per  minute. 
While  this  rate  of  drift  makes  the  INS  unsuitable  for  a  PLS,  it  does  give  an  accurate 
enough  solution  for  an  aircraft  to  navigate  safely  away  from  the  landing  area  when 
a  GPS  failure  is  detected.  This  work  will  consider  integration  of  low-,  medium-,  and 
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high-quality  INSs  as  discussed  in  Section  1.4.  The  low-quality  INS  is  used  here  to 
illustrate  the  drift  problem. 

The  accuracy  potential  of  GPS  as  a  primary  sensor  in  a  precision  landing 
system  is  evident.  However,  possible  interference  of  the  GPS  signal,  as  well  as 
possible  system  failures  from  other  sources,  remain  a  major  concern.  Whatever  the 
source  of  the  interference,  the  navigation  solution  provided  by  any  GPS  receiver 
under  interference  conditions  would  be  unreliable.  The  dependability  of  the  GPS 
signal  as  a  safety-of-flight-critical  sensor  in  a  precision  landing  system  remains  to  be 
adequately  addressed. 

1.2  Problem  Definition 

The  purpose  of  this  thesis  is  to  develop  an  effective  method  of  receiver  au¬ 
tonomous  integrity  monitoring  (RAIM)  for  GPS-aided  INS-based  PLS.  The  RAIM 
scheme  to  be  developed  in  this  work  will  accurately  indicate  the  integrity  of  the 
GPS  (and  overall  system)  navigation  solution  during  interference/jamming  or  spoof¬ 
ing  conditions. 

1.3  Summary  of  Current  Knowledge 

Recent  research  in  the  area  of  GPS-based  PLS  conducted  by  the  FAA  and 
other  researchers  has  focused  on  developing  and  verifying  a  GPS-based  navigation 
architecture  [3,10,33,37,46]  which  can  provide  the  necessary  navigational  accuracy 
for  precision  landings.  The  work  of  Gray  [10]  (using  standard  GPS)  and  Britton  [3] 
(using  DGPS)  at  AFIT  has  shown  that  an  integrated  GPS-aided  INS-based  PLS 
meets  FAA  requirements  for  Category  I  and  II  precision  approaches.  Paielli,  Russell, 
and  others  [33]  have  demonstrated  the  increased  accuracy  potential  which  may  be 
obtained  by  using  carrier  phase  GPS  measurements. 

Much  research  effort  has  focused  on  developing  a  working  RAIM  method  for 
GPS.  Comparatively  little  work  has  focused  on  integrity  monitoring  for  a  GPS-aided 
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INS-based  navigation  system.  Vasquez  [41,42]  has  applied  generalized  likelihood 
ratio  and  chi-square  testing  schemes  to  the  problem  of  interference /jamming  and 
spoofing  detection  for  the  GPS  signal  in  a  GPS-aided  INS-based  system. 

1.4  Assumptions 

The  integrated  system  to  be  used  will  include  a  ground-based  pseudolite  which 
will  provide  GPS-like  transmissions  from  a  surveyed  location  near  the  runway  [10]. 
The  main  advantage  of  a  ground-based  pseudolite  is  the  improvement  it  makes  in 
the  geometry  of  the  GPS  problem.  GPS  satellites  are  located  strictly  above  the 
user’s  horizon.  The  pseudolite  provides  a  range  measurement  from  beneath  (below 
the  horizon)  the  user.  In  this  research,  it  will  be  assumed  that  differential  GPS, 
pseudolite,  radar  altimeter,  baro  altimeter,  and  INS  navigation  signals  will  all  be 
available.  The  benefits  (to  failure  identification)  of  the  pseudolite  and  radar  altimeter 
will  be  analyzed  by  conducting  some  performance  simulations  without  one  or  both 
of  these  measurements.  Three  different  accuracies  of  INSs  will  be  used,  a  high- 
quality  INS  as  might  be  used  by  the  military,  a  medium-quality  INS  like  those 
used  in  commercial  transportation,  and  a  low-quality  INS  which  may  be  inexpensive 
enough  to  be  used  in  small  private  aircraft.  Radar  altimeter  measurements  are 
assumed  to  be  reflected  from  the  WGS-84  reference  ellipsoid  (no  terrain  effects  will 
be  modeled)  when  the  aircraft  is  at  less  than  3000  feet  AGL  (above  ground  level). 
A  radar  altimeter  is  considered  to  be  a  reasonable  part  of  such  a  PLS  because  of 
its  common  availability  on  military  aircraft  and  commercial  airliners,  as  well  as  the 
non-prohibitive  cost  of  incorporating  it  into  small  civil  aviation  vehicles. 

1.5  Literature  Review 

This  section  reviews  the  current  literature  relative  to  the  GPS-aided  INS-based 
Precision  Landing  System  (PLS).  It  also  reviews  the  literature  pertaining  to  several 
of  the  most  viable  techniques  of  failure  detection  and  identification.  The  algorithm 
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discussions  presented  in  this  section  are  general  and  quite  terse;  a  detailed  technical 
description  of  each  algorithm  is  taken  up  in  Chapter  2. 

1.5.1  Basis  of  a  GPS-aided  INS-based  PLS.  Much  recent  research  has 
focused  on  developing  and  flight  testing  a  GPS-based  PLS  to  replace  the  current 
ILS  (and  recently  cancelled  MLS)  system  [3,10,12,33].  To  ensure  safety  of  flight, 
aircraft  require  an  internal,  independent  navigation  reference  system  that  may  be 
used  when  the  GPS  signal  is  unhealthy  or  jammed.  Where  the  cost  is  justified, 
inertial  navigation  systems  (INSs)  are  used  to  provide  this  independent,  internal 
navigation  reference.  Military  and  large  commercial  aircraft  use  INSs.  The  cost 
of  lower-quality  INSs  (and  of  GPS  systems)  has  become  realistic  for  civil  aviation 
applications.  Additional  navigation  devices  may  be  required  to  aid  the  INS  during 
lengthy  GPS  outages. 

INSs  provide  very  accurate  (relative)  navigation  over  short  periods  of  time  but 
drift  substantially  (accumulate  position  bias)  over  time.  The  GPS  provides  accu¬ 
rate  navigation  information  over  long  periods  (no  drifting  trend)  but  can  easily  be 
lost  for  short  times  during  periods  of  highly  dynamic  aircraft  maneuvering  or  GPS 
signal  interference.  Integration  of  INS  and  GPS  navigation  data  with  other  navi¬ 
gation  aids  provides  increased  navigation  accuracy  during  benign  flight  and  allows 
reliable  navigation  in  the  face  of  aircraft  dynamics  or  intermittent  GPS  interfer¬ 
ence/jamming  [10,12,36].  The  standard  navigation  system  on  military  aircraft  is 
based  on  an  INS  aided  by  GPS  and  other  navigation  aids.  A  GPS  user  can  only 
receive  signals  from  satellites  located  above  the  horizon  (ideally,  transmitters  would 
be  available  from  below  the  horizon  as  well,  as  is  the  case  when  pseudolites  are 
present  in  the  local  area);  there  is  no  sensitivity  to  the  cardinal  direction  (azimuth) 
to  the  satellite.  This  geometry  makes  the  GPS  navigation  solution  least  precise  in 
the  vertical  direction.  To  help  overcome  this  weakness,  altitude  sensors  are  standard 
navigation  aids  for  systems  using  GPS.  An  integrated  system  consisting  of  (at  least) 
GPS,  an  INS,  and  barometric  altimeter  is  common  on  military  aircraft.  A  radar 
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altimeter  is  also  generally  present  on  military  and  commercial  aircraft  although  it  is 
not  typically  a  part  of  the  integrated  navigation  system.  For  this  work  a  radar  al¬ 
timeter,  when  present,  will  be  modelled  as  a  part  of  the  integrated  navigation  system. 
The  integration  of  the  radar  altimeter  was  first  done  by  Gray  [10]  and  Britton  [3] 
to  improve  the  vertical-channel  navigational  precision  during  approach  and  landing 
maneuvers. 

Where  an  integrated  navigation  system  is  in  place,  especially  in  military  appli¬ 
cations,  it  is  desired  to  develop  a  PLS  based  on  the  entire  available  set  of  navigation 
components.  Past  AFIT  theses  have  focused  on  the  development  of  a  GPS-aided 
INS-based  landing  system.  This  research  has  shown  that  a  differential  GPS-aided 
INS-based  PLS  provides  Category  I  and  II  landing  precision  [3,10].  The  reliabil¬ 
ity  of  this  landing  configuration,  especially  in  possibly  jammed  or  spoofed  areas  of 
operation,  remains  to  be  fully  investigated. 

1.5.2  Integrity  Monitoring.  This  work  focuses  on  the  problem  of  detect¬ 
ing  radio  frequency  (rf)  interference  which  would  corrupt  local  GPS  signals  and  the 
resulting  navigation  solution.  Civilian  GPS  users  will  most  often  face  benign  in¬ 
terference  from  rf  sources  such  as  microwave  or  television  transmitters.  Military 
GPS  users  must  anticipate  malignant  interference/] amming  or  spoofing  from  enemy 
sources.  The  signal  fault  identification  process  is  referred  to  as  Fault  Detection  and 
Identification  (FDI).  This  section  reviews  the  literature  pertaining  to  four  techniques 
of  FDI,  which  have  been  or  could  be  applied  to  the  current  integrity  monitoring  prob¬ 
lem.  Kalman  filtering  will  be  used  in  all  of  the  cases  examined  to  provide  integration 
of  the  INS  with  other  systems.  The  measurement  residual  signals  generated  by  the 
Kalman  filters  will  be  used  to  perform  the  FDI  functions  [13-15,20,25].  See  May- 
beck  [21-23]  for  a  thorough  presentation  of  Kalman  filter  theory  and  applications. 

1.5.2. 1  Integrated  Navigation  and  FDI  Concepts.  Real-world  naviga¬ 
tion  sensors  cannot  produce  an  exact  navigation  solution;  they  give  a  measurement. 
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corrupted  by  noise,  of  the  unknown  navigation  parameters.  Modern  integrated  nav¬ 
igation  is  implemented  by  using  Kalman  filters  to  combine  measurements  optimally 
from  multiple  information  sources,  based  upon  the  known  (or  determined)  relative 
precision  of  each  source.  Conceptually,  Kalman  filters  estimate  the  system  state, 
propagate  that  estimate  forward  to  the  next  measurement  sample  time  using  an 
assumed  system  model,  then  update  that  state  estimate  using  actual  measurement 
information.  The  difference,  at  each  update  time,  between  the  filter’s  state  predic¬ 
tion  (based  on  propagating  its  system  model)  and  the  measurement  actually  taken, 
is  called  the  measurement  residual.  The  measurement  residual  is  a  reliable  indicator 
of  how  well  the  filter’s  assumed  system  model  matches  the  actual  system.  When 
the  filter  model  disagrees  with  the  real  world,  the  characteristics  of  the  residuals  can 
provide  information  about  how  real  system  differs  from  the  filter’s  model. 

The  fundamental  objective  of  FDI  is  to  examine  the  available  information  in 
such  a  way  that  system  failures  can  be  detected  and  identified.  Navigation  sensor 
failures  are  identifiable  as  discrepancies  in  the  solutions  provided  by  different  infor¬ 
mation  sources  [44],  either  as  directly  viewed  from  the  sensor,  or  as  observed  in  the 
character  of  the  measurement  residual  of  the  integration  Kalman  filter.  The  PLS 
under  examination  here,  consisting  of  (at  least)  GPS,  an  INS,  barometric  altimeter, 
and  radar  altimeter,  is  well-suited  to  FDI  techniques  based  on  multiple,  redundant 
information  sources.  The  solution  to  the  current  problem  will  involve  checking  the 
received  GPS  navigation  data  against  the  internal  INS  navigation  solution,  which  is 
not  subject  to  external  interference.  It  is  desired  to  incorporate  the  received  GPS 
signals  when  they  are  currently  valid,  because  the  GPS  solution  accuracy  greatly 
contributes  to  the  precision  required  for  ILS-like  aiding. 

Each  FDI  algorithm  examined  differs  in  function  and  complexity.  Of  the  al¬ 
gorithms  discussed,  chi-square  testing  only  detects  the  existence  of  failures  in  the 
system,  while  the  other  three  FDI  methods  perform  some  degree  of  fault  identifica¬ 
tion  and  recovery  [22,26,40,44]. 
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1.5. 2. 2  Redundancy  (The  Voting  Method).  Perhaps  the  simplest 
and  most  reliable  failure  detection  technique  is  the  use  of  duplicate  hardware  for 
voting  [41,42,44].  Simple  voting  requires  three  redundant  information  sources.  The 
outputs  of  each  identical  element  are  compared,  allowing  them  to  vote  on  what 
information  is  valid.  If  two  elements  agree  on  the  aircraft  position  and  the  third 
provides  a  different  solution,  the  deviant  element  is  declared  to  have  failed. 

The  voting  method  for  triply  redundant  devices  has  two  restrictive  weaknesses. 
Once  any  single  element  has  failed  and  is  removed  from  the  redundant  triad,  the  sys¬ 
tem  can  detect,  but  is  unable  to  isolate,  a  second  failure.  If  a  second  failure  occurs, 
the  system  cannot  determine  which  element  is  in  error  and  both  remaining  elements 
should  be  taken  off  line  in  order  to  maintain  the  integrity  of  the  system’s  naviga¬ 
tion  solution.  A  second  weakness  of  the  voting  method  is  the  additional  expense, 
space,  weight,  and  computational  ability  required  for  redundant  hardware  [26,40,44]. 
The  requirement  for  these  additional  resources  typically  makes  the  voting  method 
unacceptable  for  avionics  applications. 

1.5. 2. 3  Chi-Square  Testing.  The  chi-square  testing  algorithm  [40, 
44]  assumes  the  use  of  a  Kalman  filter  to  combine  information  from  the  available 
navigation  sensors.  This  method  is  based  on  monitoring  the  measurement  residuals  of 
the  Kalman  filter.  The  filter’s  residuals  and  their  filter-computed  residual  covariance 
are  monitored.  Measurement  residuals  for  a  filter  having  the  correct  model  of  the  real 
world  should  display  four  well-defined  characteristics.  Residuals  should  be  white, 
Gaussian,  zero  mean,  and  have  covariance  HP~ -|-  R  (or  closely  approximate 
those  properties)  [6,13,21,22,40].  Gaussian-ness,  zero-mean- ness,  and  covariance 
HP~  +  R  are  exploited  in  the  test  conducted  in  this  fault  detection  method.  The 
chi-square  test  declares  a  failure  when  the  tested  properties  of  the  filter  residuals  are 
inconsistent  with  those  expected  from  a  filter  having  the  correct  model.  For  instance, 
one  can  run  a  hypothesis  test  (to  some  confidence  level)  that  95%  of  a  residual’s 
scalar  component  lies  between  0±  twice  the  filter-computed  standard  deviation  for 
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that  residual  component.  The  chi-square  test  is  easy  to  implement  and  runs  quickly 
but  provides  only  a  binary  fail/no-fail  indication  of  system  operation.  It  is  not  useful 
for  failure  identification  [41,42,44]. 

1.5. 2. 4  Generalized  Likelihood  Ratio  Testing.  Like  chi-square  testing. 
Generalized  Likelihood  Ratio  testing  (GLR)  uses  the  measurement  residuals  of  the 
navigation  system’s  Kalman  filter  [40-42,44,45].  The  GLR  algorithm  attempts  to 
detect  and  isolate  failures  by  knowing  the  effect  that  each  failure  has  on  the  character 
of  the  filter  residuals.  Maximum  Likelihood  Estimation  (MLE)  is  used  to  determine 
which  condition  is  most  probable,  given  the  current  character  of  the  filter  residuals. 
In  general,  to  implement  the  GLR  method,  the  designer  is  required  to  provide  a 
signature  failure  matrix  which  provides  the  algorithm  with  the  description  of  how 
different  failures  modify  the  filter  residuals  [44].  Uncertain  parameters  such  as  failure 
magnitudes  may  also  be  estimated  by  modelling  their  effects  in  the  signature  failure 
matrix. 

Van  Trees  [40]  shows  the  development  of  multiple  GLR  testing.  In  multiple 
GLR  testing,  the  system  Kalman  filter  is  designed  based  on  the  no-fail  condition. 
Matching  filters  are  designed  based  on  the  expected  failure  conditions  (including 
no-fail).  Each  matching  filter  generates  a  GLR  using  MLE.  Based  on  the  GLRs, 
central  testing  logic  determines  which  matching  filter  best  matches  the  real  world. 
Once  a  failure  is  detected,  a  corrective  signal  can  be  fed  back  to  the  Kalman  filter 
for  adaptation  to  the  failure.  Further  detail  on  GLR’s  is  presented  by  Willsky  and 
Jones  [44,45].  See  Vasquez  [41,42]  for  an  example  of  multiple  GLR  testing  applied 
to  a  navigation  problem. 

1.5. 2. 5  MMAE.  The  MMAE  technique  uses  multiple  Kalman  filters 
running  in  parallel  to  model  the  dynamics  of  the  system  (in  this  case  the  PLS)  under 
different  conditions  of  failed  or  no-fail  operation  [1,20,22].  A  separate  Kalman  filter 
is  designed  for  each  failure  condition  and,  during  operation,  the  residuals  are  used 
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to  determine  how  well  each  filter  models  the  current  system  state.  The  MMAE 
technique  is  similar  to  multiple  GLR  testing  in  many  ways,  but  differs  in  its  structure 
and  decision  making  methods.  Multiple  GLR  testing  uses  a  variety  of  matching 
filters  to  model  failure  characteristics  of  the  residuals  of  a  single  Kalman  filter.  In 
the  MMAE  algorithm,  a  probability  of  model  accuracy  ranging  from  zero  to  one 
is  computed  for  each  elemental  filter  within  the  MMAE  structure.  Each  separate 
filter’s  state  estimates  are  scaled  by  the  computed  probabilities  that  each  filter  has 
the  correct  model  of  the  real  world.  These  scaled  state  estimates  are  then  added 
together  to  produce  a  probability-weighted  blending  of  the  state  estimates.  This 
blending  technique  has  particular  advantages  when  the  real  world  system  is  not 
exactly  described  by  any  one  filter  model,  but  exists  in  the  parameter  space  between 
discrete  filter  models.  In  such  a  situation  the  filter  residuals  of  more  than  one  of  the 
modeled  conditions  look  reasonably  valid.  The  weighted  state  estimates  provided  by 
the  MMAE  algorithm  are  used  to  construct  the  navigation  solution  and  to  provide 
system  feedback.  A  major  strength  of  the  MMAE  technique  lies  in  its  capacity  for 
rapid  and  valid  adaptation.  Multiple  filters  are  running  in  parallel,  making  multiple 
sets  of  residuals  available  at  all  times  with  which  to  decide  which  filter  model  looks 
best  (according  to  the  “good”  residual  qualities  [21,22]  described  in  Section  1.5. 2. 3). 

For  robust  operation,  the  bank  of  Kalman  filters  must  be  specified  in  such 
a  way  that  it  spans  the  entire  failure  space  of  the  application  of  interest.  Rigor¬ 
ous  general  convergence  proofs  for  MMAE  do  not  exist  [7,9,43],  but  some  limited 
proofs  have  been  established  [4,5,22,24]  and  experience  has  shown  that,  as  long  as 
the  failure  condition  of  the  system  remains  within  the  span  of  the  filter  bank,  the 
MMAE  technique  is  robust  to  that  condition  and  will  respond  to  it  very  quickly  and 
accurately. 

1.5.3  Literature  Review  Conclusion.  If  it  is  to  be  used  in  actual  opera¬ 
tion,  the  proposed  GPS-aided  INS-based  PLS  will  require  rapid  and  accurate  fault 
detection  and  identification  to  establish  the  integrity  of  the  PLS  aid  at  all  times. 
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The  error  detection  techniques  presented  in  Section  1.5.2  are  representative  of  the 
methods  available  to  accomplish  the  required  fault  detection. 

The  redundancy  technique  is  not  feasible  for  this  application.  Multiple  GPS 
receivers  on  a  single  aircraft  would  probably  not  be  restrictive  on  the  basis  of  size, 
weight,  or  computational  loading,  but  each  receiver  would  be  subject  to  the  same 
failure  environment,  e.g.  interference/jamming  or  spoofing,  and  so  nothing  would 
be  gained  by  this  redundancy. 

Because  of  its  speed  and  simplicity,  chi-square  testing  has  been  widely  used  to 
detect  failures.  However,  chi-square  testing  alone  has  no  ability  to  isolate  failures 
and  so  does  not  have  much  applicability  to  the  integrity  monitoring  problem  under 
examination,  if  used  alone.  This  testing  method  has  been  effectively  used  when 
combined  with  another  test  to  perform  the  identification  function. 

Multiple  GLR  testing  and  MMAE  are  both  well  suited  to  the  FDI  problem. 
Both  of  these  techniques  are  flexible  enough  to  be  applied  to  the  complex  GPS-aided 
INS-based  problem.  This  work  uses  the  MMAE  technique  to  perform  the  failure 
detection  and  identification.  Vasquez  [41,42]  has  used  a  combination  of  chi-square 
and  GLR  testing  to  detect  and  isolate  failures  in  a  GPS-aided  INS-based  navigation 
system.  Performance  comparisons  will  be  made,  as  appropriate,  between  the  MMAE 
FDI  and  GLR/chi-square  FDI  algorithms. 

1.6  Scope 

This  work  will  be  limited  to  the  development  of  a  multiple  model  adaptive 
estimation  (MMAE)  architecture  to  provide  the  integrity  monitoring  previously  de¬ 
scribed.  Comparison  will  also  be  made  to  the  generalized  likelihood  ratio  and  chi- 
square  testing  algorithms  developed  by  Vasquez.  The  failure  modes  to  be  addressed 
will  be  limited  to  the  onset  of  interference/jamming  and  spoofing.  Filters  based 
on  ramp  failures  will  not  be  included  in  the  working  MMAE  filter  bank.  When  a 
ramp  occurs,  the  elemental  filter  with  the  hypothesized  constant  parameter  value 
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that  most  closely  matches  the  current  ramp  value  will  receive  the  highest  probabil¬ 
ity  weighting.  The  ramp  should  be  detectable  as  a  growing  trend  in  the  MMAE’s 
blended  estimate  of  the  affected  parameter  (and  a  gradual  shift  in  the  computed 
probabilities  of  each  hypothesized  parameter  being  correct). 

1.7  Methodology  Overview 

The  research  consisted  first  of  studying  the  stochastic  and  dynamics  modeling 
of  the  GPS-aided  INS-based  navigation  system,  primarily  as  represented  in  the  theses 
of  Gray  [10],  Britton  [3],  and  Negast  [30].  The  model  development  work  accomplished 
in  these  previous  AFIT  theses  was  duplicated  and  then  confirmed.  Once  the  model 
had  been  successfully  duplicated,  the  research  effort  focused  on  failure  detection 
within  the  augmented  system.  The  failure  detection  analysis  of  the  multiple  model 
architecture  used  was  accomplished  using  Multiple  Model  Simulation  for  Optimal 
Filter  Evaluation  (MMSOFE  [32]).  MMSOFE  allows  for  the  simultaneous  testing  of 
multiple  extended  (or  linear)  Kalman  filter  (EKF)  models  in  an  MMAE  architecture. 
MMSOFE  is  based  on  the  MSOFE  [28]  program  designed  for  the  testing  of  single 
Kalman  filters.  MSOFE  was  written  in  FORTRAN  77  [39]  and  its  use  involves 
significant  modification  of  up  to  14  user-definable  modules. 

Comparison  to  the  work  of  Vasquez  is  performed.  Vasquez  used  general¬ 
ized  likelihood  ratio  (GLR)  and  chi-square  testing  schemes  to  detect  and  estimate 
interference/! amming-  spoofing-induced  failures  in  a  GPS-aided  INS-based  sys¬ 

tem.  GLR  techniques  are  based  on  residual  monitoring  and  can  be  very  effective. 
The  failure  detection  and  identification  results  of  Vasquez’  work  are  compared  to 
those  generated  using  the  MMAE-based  FDI  methods. 
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1.8  Overview  of  Thesis 

Chapter  2  presents  the  theory  used  in  this  research.  Kalman  filter  theory  is 
introduced,  along  with  the  basics  of  several  FDI  algorithms,  including  MMAE,  GLR, 
and  chi-square  testing. 

Chapter  3  describes  the  navigation  system’s  parameters  and  operation  through 
an  overall  system  description.  Models  for  each  of  the  PLS  components,  including 
the  INS  (with  barometric  altimeter),  GPS,  DGPS,  pseudolite,  and  radar  altimeter, 
are  defined  in  detail.  The  failure  models  used  to  represent  interference /jamming 
and  spoofing  failures  are  developed,  along  with  a  description  of  the  MMAE-based 
methods  that  are  used  to  detect  those  failures. 

Results  of  the  work  done  are  shown  in  Chapter  4,  including  an  analysis  of  the 
FDI  (and  navigation  guidance)  performance  observed.  Chapter  5  summarizes  the 
research  through  conclusions  and  recommendations. 
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2.  Theory 


2.1  Overview 

This  chapter  presents  the  fundamental  theory  of  the  Kalman  filter,  and  the 
sampled-data  Kalman  filter  equations  used  in  this  work  are  developed.  The  ex¬ 
tended  Kalman  filter  (EKF),  a  well-established  ad  hoc  method  based  on  the  Kalman 
filter,  is  used  in  estimation  problems  involving  nonlinear  dynamics  and/or  nonlinear 
measurement  models.  The  equations  defining  (EKF)  are  presented  in  this  chapter, 
along  with  a  more  detailed  discussion  of  the  MMAE,  GLR,  and  chi-square  failure 
detection  methods  introduced  in  Chapter  1. 

2.2  The  Extended  Kalman  Filter 

The  EKF  allows  for  nonlinear,  time- varying  system  dynamics  and/or  measure¬ 
ment  vectors,  as  are  found  in  this  GPS/INS  navigation  problem.  In  simple  linearized 
Kalman  filtering  (LKF),  the  dynamics  and/or  measurement  equations  are  linearized 
through  first-order  perturbation  techniques  about  a  fixed  nominal  trajectory.  The 
LKF  is  the  conceptual  basis  for  the  EKF.  During  operation  the  EKF  is  continually 
relinearized  about  the  most  recent  state  estimate  trajectory  rather  than  about  a 
fixed  nominal  trajectory. 

For  the  sampled  data  Kalman  filter,  let  the  system  model  be  expressed  as  a 
state  equation  of  the  form 


x(f)  =  f[x(t),t]  +  G{t)w(f)  (2.1) 

where  the  state  dynamics  vector  f[x(f),f]  is  a  (possibly)  nonlinear  function  of  the 
state  vector  x(i)  and  time  t.  Let  w(i)  be  a  white  Gaussian  noise  with  mean 

E[w{t)]  =  0  (2.2) 
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and  noise  strength  Q(t)  defined  by 


E[wit)w^{t  +  T)]  =  Qit)6{T)  (2.3) 

The  discrete-time  measurements,  z(ti),  are  modeled  as  a  (possibly)  nonlinear 
vector  of  functions  of  the  state  vector  and  time,  h[x(ti), i,],  plus  additive  white 
Gaussian  noise: 

z{ti)  =  h[x(ti),fi] -h  v(ti)  (2.4) 

where  h[x(tj),  tj]  is  the  nonlinear  observation  vector  and  v(ti)  is  a  zero- mean  discrete¬ 
time  white  Gaussian  noise,  independent  of  the  dynamics  driving  noise  w(t)  and 
having  covariance  R(ti)  defined  by 

R(ti)  for  U  =  tj 
0  for  ti  tj 

The  LKF  is  based  on  perturbation  states  about  a  nominal  state  trajectory  x„(t) 
satisfying  x„(to)  =  and 

x„(t)  =  f[x„(t),t]  (2.6) 

using  the  same  f  [•,  •]  as  in  Equation  (2.1).  The  nominal,  noise-free  measurements  are 
also  based  on  the  nominal  states  and  are  defined  as 

z„(ti)  =  h[x„(ti),ti]  (2.7) 

The  perturbation  states  are  found  by  subtracting  the  nominal  states  in  Equation  (2.6) 
from  the  original  states  in  Equation  (2.1): 

[x{t)-Xn{t)]  =  i[x{t),t]-f[xn{t),t]  +  G{t)w{t)  (2.8) 


(2.5) 


F;[v(ti)v^(ti)]  = 
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Equation  (2.8)  is  approximated  to  first  order  through  a  truncated  Taylor  series  ex¬ 
pansion  (letting  (5x(t)  be  a  first-order  approximation  to  [x(t)  Xn(t)]): 

Sx(t)  =  F[t-x{t)]Sxn{t)  +  G{t)w{t)  (2.9) 

where  Sx{t)  are  the  perturbation  states.  The  definitions  for  G{t)  and  w{t)  are 
unchanged,  and  the  new  linearized  dynamics  matrix  F[t;x„(t)]  is  found  by  taking 
partial  derivatives  of  f[x(t),t]  with  respect  to  x(t)  and  evaluated  at  the  nominal 
values  for  the  trajectory  x„(t): 

F[<;x„(t)]  =  5^1^  (2.10) 

The  discrete-time  perturbation  measurements  are  similarly  approximated  to  first 
order  from  the  measurement  difference  equation 


6z{ti)  =  z{ti)-Zniti)  =  h[x{ti),ti]-h[Xn{U),ti]+v{ti)  (2.11) 


yielding  the  perturbation  form  (letting  6z{ti)  be  a  first-order  approximation  to 

[z(tj)  Zn(ti)]. 

Sz{ti)  =  ll[ti-,x{ti)]6xn{ti) +  (2.12) 

The  same  partial  derivative  methods  used  to  derive  the  linearized  state  dynamics 
matrix  F[t;x„(t)]  are  used  again  to  derive  the  linearized  observation  matrix: 


X7i(t j)] 


5h[x(fi),ti] 


dx 


X 


(2.13) 


Because  the  LKF  generates  error  state  estimates,  8x{t),  they  must  be  added  to  the 
nominal  states  to  provide  whole  state  estimates  x{t)  in  the  form 


x{t)  =  Xn{t)  +  6x{t) 


(2.14) 
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If  the  nominal  and  the  “true”  trajectories  differ  too  greatly,  a  linearized  Kalman 
filter  will  yield  erroneous  state  estimates  because  the  condition  for  neglecting  higher 
order  terms  of  the  Taylor  series  expansion  is  violated.  The  EKF  will  be  formed  by 
linearizing  about  the  most  recent  state  estimate  x  rather  than  about  the  nominal 
trajectory  x„,  as  is  done  in  the  LKF.  The  following  sampled  data  FKF  equations  use 
the  notation  t\ti  to  represent  the  value  of  a  given  variable  at  time  t,  conditioned  on  the 
measurements  taken  through  time  U.  Also,  ti~  represents  the  value  after  propagation 
but  prior  to  the  measurement  update  at  sample  time  t,-,  and  corresponds  to  the 
value  after  the  measurement  update.  The  state  estimates  x(t|ti)  and  covariance 
values  are  propagated  from  ti  to  U+i  by  solving  the  following  differential 

equations: 


St{t\ti)  =  f[x(t|ii). 

i] 

(2.15) 

=  F[t;  x(t  |fO]P(i  mu)]  +  G(t)Q(t)G^(t) 

(2.16) 

where 

F[,xW,.)l  = 

(2,17) 

X  =  x(t\u) 

and  initial  conditions  are  given  by: 

m\u)  =  x(^i+) 

(2.18) 

P{U\U)  =  P(ti+) 

(2.19) 

The  discrete-time  measurements  are  processed  in  the  FKF  through  the  update  equa¬ 
tions: 

K(ti)  =  P(tr)H^[ii;x(tr)]{H[ti;x(tr)]P(^r)H^[ii;x(ir)] +  R(ti)}  ^  (2.20) 


x(i,.+)  = 

m  )  +  K(ti)  {zi  -  h[x(ti 

(2.21) 

P(l,.+)  = 

P(fr)-K(ti)H[i,;x(fr)]P(t-) 

(2.22) 
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where 


P(ii  )  =  P(iilii_i)  and  x(ii  )  =  (2.23) 


as  produced  by  the  most  recent  propagation  cycle.  The  variable  Zj-  represents  the 
actual  realization  of  the  measurement  z(ti),  )]  is  given  by: 


dh[x{ti),ti 


dx 


X 


(2.24) 


and  K(ti)  is  the  discrete-time  Kalman  filter  gain.  Recall  from  the  previous  page  that, 
for  the  EKF,  the  measurement  and  dynamics  vectors  are  calculated  about  the  last 
state  estimate  x{ti~)  and  the  state  trajectory  x(t|t8),  rather  than  about  the  nominal. 

Before  proceeding  to  a  discussion  of  failure  detection  and  identification  (FDI) 
techniques,  the  reader  should  make  note  of  Equation  (2.21).  The  term  {z^  —  h[x(f,“),  fi]} 
on  the  right  hand  side  of  Equation  (2.21)  is  called  the  measurement  residual.  At 
each  sample  time  fj,  the  residual  is  the  difference  between  the  actual  measurement  of 
the  real  world  z(t,)  and  the  filter’s  best  prediction  of  what  the  measurement  should 
have  been  before  it  arrived,  h[x(fj“),fi].  The  filter’s  prediction  is  based  on  its  model 
of  the  system,  so  the  characteristics  of  the  measurement  residual  are  an  indication 
of  how  well  the  filter’s  model  currently  matches  the  real  world.  The  residuals  from 
a  filter  (the  filter  of  a  number  of  possible  filters)  having  the  correct  model  will 
be  (for  a  linear  KF,  or  to  first  order  for  an  EKF)  white,  Gaussian,  zero-mean,  and 
have  covariance  [Hfc(fj)Pfc(fi“)H*;^(fj)  -t-  R*(fi)].  This  information  contained  in  the 
residual  is  the  basis  of  the  FDI  techniques  discussed  in  the  next  section. 


2.3  Failure  Detection 

The  MMAE  and  the  GLR/chi-square  methods  of  FDI  [22,  40]  are  residual- 
based  techniques  built  upon  the  Kalman  filter  just  developed  in  Section  2.2.  The 
basic  ideas  upon  which  these  FDI  methods  are  based  were  discussed  in  Section  1.5.2. 
This  section  presents  a  theoretical  discussion  of  the  MMAE,  GLR,  and  chi-square 
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algorithms.  See  VanTrees  [40]  for  a  more  complete  development  of  the  GLR  and  chi- 
square  methods.  Maybeck  [22]  presents  a  clear  and  rigorous  development  of  MMAE 
theory. 


2.3.1  Multiple  Model  Adaptive  Estimation.  Multiple  Model  Adaptive  Esti¬ 
mation  (MMAE)  can  provide  simultaneous  state  and  parameter  estimation.  MMAE 
is  composed  of  multiple  “elemental”  Kalman  filters  running  in  parallel,  each  using 
an  identical  deterministic  input  and  measurement  environment.  Each  of  the  individ¬ 
ual  elemental  filters  may  model  a  different  set  of  system  parameter  values  based  on 
known  possible  operating  conditions,  or  may  model  some  possible  failure  condition 
such  as  greatly  increased  measurement  noise  covariance  R(fi),  which  might  represent 
the  effects  of  GPS  interference/] amming.  The  magnitude  of  a  failure  can  also  be  esti¬ 
mated  by  using  multiple  elemental  filters  with  different  assumed  failure  magnitudes. 
The  residuals  from  each  filter  are  used  to  calculate  the  conditional  probabilities  that 
each  of  the  filters  has  the  most  correct  model.  Uncertain  parameters  are  estimated 
in  this  way  because  specific  parameter  values  are  associated  with  each  filter.  The 
conditional  probabilities,  also  referred  to  as  hypothesis  conditional  probabilities,  will 
be  used  as  weights  for  blending  the  state  estimates  from  the  elemental  filters  to 
produce  the  final  blended  MMAE  state  estimates.  See  Figure  2.1  for  a  graphical 
representation  of  the  MMAE  algorithm. 

The  conditional  probability  pk{ti)  for  the  elemental  filter,  k  =  1,2,. . .  ,K 
is  determined  by: 


Pk{ti)  = 


12j=l  /z(«,)|a,Z(<,_i)(Zi|aj, 


(2.25) 


where,  for  this  development,  the  parameters  may  represent  measurement  bias 
magnitudes  and/or  measurement  noise  covariance  values.  The  numerator  of  Equa¬ 
tion  (2.25)  is  the  product  of  two  terms.  The  first,  pk{ti-i),  is  merely  the  most  recent 
value  of  the  conditional  probability  for  the  elemental  filter,  making  this  equation 
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Figure  2.1  Multiple  Model  Adaptive  Estimation 


iterative.  All  K  numerator  terms  for  the  K  elemental  filters  at  must  be  available 
before  the  denominator  and  consequently  the  conditional  probability  at  the  cur¬ 
rent  time,  ti,  can  be  calculated.  The  first  numerator  term  is  the  probability  density 
function  of  the  current  measurements,  conditioned  on  both  the  assumed  parameter 
values  and  the  observed  past  measurements.  This  probability  density  function  is 
computed  by: 


/z(<0|a,ZPi_a)(z.|a„Zi_i)  =  l/(27r)  2  | Afc(ti) l^exp  {•} 

where  m  is  the  measurement  vector  dimension,  and  the  filter  residual  is  given  (for 
linear  filters)by: 


rk{ti)  =  [z{ti)  -  Hk{ti)9.k{ti  )] 


(2.27) 
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and 

M'^i)  =  {[z(^i) -hfc[xfc(ir),  (^i)]}  (2.28) 

for  the  more  general  extended  Kalman  filters,  and  where  the  residual  covariance  is 
computed  by  the  elemental  filter  as: 

Ak{ti)  -  [Hk{ti]Fk{tniik^{U)  +  Rk{ti)]  (2.29) 

The  filter  residual,  rfc(t,),  is  dependent  upon  the  current  measurement,  z(tj),  the 
measurement  matrix,  Hfc(tj)  (or  the  measurement  function  hjt(-,  i,)  for  an  EKF),  and 
the  state  estimate  prior  to  the  measurement,  x(t,~).  The  k*^  elemental  filter’s 
state  estimation  error  covariance  matrix  before  the  measurement,  Pfc(ti~),  the 
measurement  matrix,  and  the  observation  noise  covariance  matrix,  Rfc(ti),  are  used 
to  construct  the  residual  covariance  for  each  filter.  The  Kalman  filter  equations  and 
notation  were  discussed  in  Section  2.2. 

If  the  residuals  in  the  k*^  filter  display  a  mean  of  zero  and  the  correspond¬ 
ing  filter- computed  residual  covariance  Ak{ti),  the  exponential  term  {•}  in  Equa¬ 
tion  (2.26)  is  approximately  equal  to  — — ,  where  m  is  the  dimension  of  z(fj)  and 
r(ti).  However,  if  the  residuals  are  much  bigger  than  anticipated  due  to  the  wrong 
parameter  hypothesis,  the  exponential  term  {•}  in  Equation  2.26  is  a  much  larger 
negative  number  (A^^  is  positive  definite),  so  pk  decreases  exponentially.  The  de¬ 
nominator  of  Equation  (2.25)  represents  the  sum  of  the  numerator  terms  from  each 
elemental  filter  computed  at  time  i,-.  This  scales  the  numerator  to  ensure  that  the 
sum  of  the  p^s  is  always  one.  However,  a  difficulty  arises  if  the  conditional  probabil¬ 
ity  of  a  state  estimate  were  to  become  zero.  In  this  case  the  conditional  probability, 
Pk,  becomes  identically  equal  to  zero  for  all  time  thereafter.  Once  equal  to  zero, 
the  iterative  form  of  Equation  (2.25)  locks  that  pk  at  a  zero  value  for  all  subsequent 
calculations,  even  if  that  filter  begins  to  produce  good  state  estimates.  To  avoid  this 
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problem,  a  lower  bound  (threshold)  should  be  set  on  each  pk  to  prevent  such  a  zero 
lock-in  condition  [22]. 

The  state  estimates  from  each  filter,  Xjt,  are  then  scaled  by  the  corresponding 
weighting  factor  pk-  These  weighted  state  estimates  are  summed  for  all  K  filters, 
resulting  in  the  Bayesian  blended  state  vector  estimate: 

K 

‘  (2-30) 

k=l 

The  covariance,  of  the  blended  solution  is  given  by  [22]: 

=  SPfc(ii){Pfc(ii''')  +  [xfc(ii'^)-x(fi+)][xfc(ti+)-x(fi+)]^}  (2.31) 

k=\ 

The  parameter  estimates,  a^,,  are  calculated  by  scaling  the  assumed  parameter  values 
from  each  elemental  filter  by  the  corresponding  weighting  factor,  p*,  in  the  same 
manner  as  for  the  state  estimates.  The  weighted  estimates  of  all  K  filters  are  summed 
according  to  the  relationship: 

K 

=  X!aA,-pfc(^i)  (2.32) 

k=l 

This  Bayesian  form  of  adaptive  estimation  is  depicted  in  Figure  2.1. 

As  may  be  inferred  from  the  preceding  presentation,  failure  detection  and  iden¬ 
tification  using  MMAE  becomes  virtually  automatic  once  the  estimator  has  been 
developed.  FDI  is  accomplished  by  observing  the  conditional  probabilities,  pk,  of 
the  elemental  filters.  The  filter  with  the  highest  probability  pk  is  based  on  the  model 
which  most  closely  matches  the  real  world  (truth  model  in  the  case  of  simulations). 
Most  likely  the  true  parameter  value  will  not  exactly  match  that  modeled  in  any  one 
of  the  elemental  filters,  so  the  unknown  parameters  will  be  estimated  via  blending 
in  a  manner  similar  to  the  MMAE  state  estimation. 


2-9 


This  thesis  research  is  directly  concerned  with  detecting  and  isolating  failures 
in  the  form  of  GPS  interference/] amming  or  spoofing.  Several  Kalman  filters  are 
used  which  differ  only  by  the  assumed  measurement  noise  variance  or  measurement 
signal  bias  magnitude.  Measurement  signal  biases  model  measurement  jumps  and/or 
ramps  (through  the  MMAE  blending  just  discussed)  as  might  result  from  intentional 
spoofing.  Increased  measurement  noise  variance  would  model  interference/] amming- 
type  real  world  measurement  noise.  Nominally,  for  the  case  of  measurement  signal 
biases,  two  filters  with  different  positive  bias  magnitudes,  a  matched  set  with  neg¬ 
ative  magnitudes,  and  one  filter  with  a  zero  bias  magnitude  can  be  assumed  (i.e., 
K  =  5,  where  K  is  the  number  of  elemental  Kalman  filters).  Using  variances  of  mea¬ 
surement  noise  instead  of  biases  in  another  set  of  elemental  filters  could  similarly 
model  interference/] amming  noise.  As  was  stated  in  Section  1.6,  there  is  no  need  to 
include  elemental  filters  based  on  a  ramp  failure  of  an  assumed  slope  and/or  time  of 
onset.  When  a  ramp  occurs,  the  elemental  filter  with  the  bias  value  that  most  closely 
matches  the  current  ramp  value  will  receive  the  highest  probability  weighting.  The 
ramp  should  be  observed  as  a  growing  trend  in  the  MMAE’s  blended  estimate  of  the 
measurement  bias. 

2.3.2  Moving-Bank  MMAE.  The  MMAE  algorithm  propagates  multiple 
filters  forward  in  time,  continually  selecting  the  filter,  or  weighted  combination  of 
filters,  that  has  the  best  model  of  the  real  world.  Often,  the  possible  parameter 
space  is  so  large  that  completely  discretizing  it  requires  many  more  filter  hypothe¬ 
ses  than  can  realistically  be  run  in  parallel.  In  this  case,  a  small  group  of  filters, 
whose  parameter  assumptions  are  in  the  close  neighborhood  of  the  current  parame¬ 
ter  estimate,  are  chosen  to  be  active  at  one  given  time.  If  the  estimated  parameter 
moves  significantly,  the  bank  of  active  filters  is  moved  to  be  centered  around  the  new 
parameter  estimate.  This  algorithm  is  called  a  moving-bank  MMAE. 
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2.3.3  Generalized  Likelihood  Ratio  Testing.  The  primary  goal  of  GLR  is  to 
define  a  likelihood  function  /(t^,  6)  that,  when  compared  to  a  threshold,  will  identify 
the  onset  of  a  failure  such  as  interference/jamming  or  spoofing.  The  GLR  algorithm 
is  depicted  in  Figure  2.2.  Multiple  hypotheses  are  established  with  a  Kalman  filter 


FDI  Block 


Figure  2.2  Multiple  GLR  Testing 


based  on  hypothesis  Ho  (no  failure)  and  matching  filters  based  on  hypothesis  Hk  (the 
failure  has  occurred).  The  matching  filters  do  not  provide  state  estimation  but 
are  designed  for  failure  detection  by  assuming  failures  in  the  system  to  be  modelled 
as  some  variation  in  the  actual  measurement  beyond  those  variations  caused  by  the 
dynamics  of  the  system.  Each  matching  filter  is  designed  to  inject  the  kind  of  residual 
modification  that  would  actually  be  experienced  in  the  system  Kalman  filter  if  that 
failure  (modelled  by  Hk)  had  actually  occurred.  The  failure  vector  d(tj)  is  m— by— 1 
where  m  is  the  number  of  measurements.  Non-failed  elements  are  indicated  in  d(tj) 
as  zeros,  while  I’s  in  the  failure  vector  correspond  to  a  measurement  device  assumed 
to  be  induced  with  a  failure.  The  arguments  of  the  likelihood  function  introduced 
above  are  ti  and  0.  As  may  be  seen  in  Figure  2.3,  the  variable  t,-  represents  the  time 
of  failure  onset,  while  to  represents  the  current  time,  and  to-N  the  time  of  the  trailing 
edge,  termed  the  “beginning”,  of  the  GLR  search  window.  As  the  search  window 
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Failure 


Search  Window 


Figure  2.3  GLR  Search  Window 


proceeds  forward  in  time,  the  failure  function  n{ti,  6)  is  either  1  or  0  depending  on 
whether  the  detected  failure  is  currently  before  or  after  the  9  index  value  in  the  GLR 
algorithm’s  search  window.  In  general,  6  takes  on  each  index  value  within  the  search 
window  at  each  time  step,  so  that  the  failure  may  be  detected  at  the  earliest  possible 
time  after  its  onset.  In  order  to  simplify  the  GLR  algorithm,  9  is  permanently  set  to 
N,  so  the  onset  of  a  failure  will  not  be  detected  until  it  reaches  the  beginning  of  the 
search  window.  Riggins  [35]  showed  that  this  simplification  gives  up  comparatively 
little  identification  performance  while  significantly  reducing  the  computational  load 
required  for  implementation. 

Under  nominal  conditions,  the  Kalman  filter  residuals  r°(L)  are: 

r°{ti)  =  Zi  -  h[x(L"),L]  (2.33) 

and  the  residuals  can  be  expressed  for  each  hypothesis  (matching  filter)  as 

TCo  :  r(L)  =  r°(L)  Hk  :  r(L)  =  +  m{ti,9)v  (2.34) 

For  a  Kalman  filter  successfully  tracking  the  true  states,  the  nominal  no-failure  resid¬ 
ual  r°(L)  will  appear  as  zero-mean  white  Gaussian  noise  of  covariance 
[Hfc(L)Pfc(ti“)Hfc^(ii)  +  Rfc(L)]-  When  a  failure  is  induced  on  the  measurements,  a 
signal  m.{ti,9)i>  will  also  be  present  in  the  residuals,  where  u  is  the  unknown  mag¬ 
nitude  and  m(ti,  9)  will  be  presented  momentarily.  The  goal  of  the  GLR  algorithm 
is  to  identify  this  signal  by  recognizing  variations  in  the  residuals  from  their  nor- 
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mal  unfailed  values.  The  GLR  tests  are  particularly  good  at  detecting  jumps  in  the 
residuals,  with  the  key  being  how  closely  the  matching  filters  model  actual  failures. 
The  GLR  algorithm  is  a  function  of  overall  system  dynamics  and  behavior  ($,  H, 
and  Kalman  filter  gain  K);  this  is  shown  mathematically  below.  The  failure  residual 
offset  m(ti,  0)  is  found  through 

xn{%e)  =  H(ti)y(ti,6>)  +  (2.35) 

where  the  recursive  failure  quantity  y(ijq.i,^)  is  given  by 

y(L+i,^)  =  ^{ti+i,ti)[l-K{ti)ll{ti)]y{ti,0)-^{ti+i,ti)K{ti)d{ti)n{ti,9)  (2.36) 

and  $(ti+i,tj)  is  the  state  transition  matrix  for  the  dynamic  system  model’s  prop¬ 
agation  from  sample  time  U  to  time  U+i.  If  the  failure  is  assumed  to  occur  at  the 
beginning  of  the  GLR  search  window  (i.e.  n(tj,  ^)  =  1  for  all  t,),  the  above  equations 
can  be  simplified  to 

m{U)  =  H(ti)y(L)  +  d{ti)  (2.37) 

y{ti+i,0)  =  ^{ti+i,ti)[I-K{ti)li{ti)]y{ti,e)  -  ^{ti+i,ti)K{ti)d{ti)  (2.38) 

The  consequence  of  this  simplification  is  a  delay  in  detecting  a  failure  because  it  is  not 
realized  until  it  reaches  the  beginning  of  the  GLR  search  window.  The  combination  of 
the  Kalman  filter  outputs  and  the  matching  filter  model  will  determine  the  magnitude 
of  the  likelihood  function  defined  as; 


=  s^{ti,e)c-\ti,e)s{ti,e)  (2.39) 

where 

S{tuO)  =  ^m^(G,6')A-^(G)r(G)  (2.40) 
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given 


C{ti,e)  =  (2.41) 

j  =  l 

Afe)  =  H((i)P(;r)H(i,)  +  R(i,.)  (2.42) 


and  the  maximum  likelihood  estimate  of  the  unknown  magnitude  of  the  failure,  z/, 
is  found  by 

The  residual  covariance  A(tj)  and  the  residuals  are  combined  with  Equations  (2.35) 
and  (2.36)  or  Equations  (2.37)  and  (2.38)  to  give  a  linear  combination  of  the  residuals 
S(ti,0)  and  a  deterministic  value  C{ti^6)  defined  in  Equations  (2.40)  and  (2.41) 
above.  Finally,  a  decision  rule  based  on  a  threshold,  e,  would  be 


l{ti,0)  >  e  =>  Declare  FAILURE 
l(ti,0)  <  e  =4>  Declare  NO  FAILURE 


(2.44) 


2.3.4  Chi-Square  Testing.  A  chi-square  test  is  based  on  the  Kalman  filter 
residuals  v{tj)  which  are  zero-mean  and  white  with  known  residual  covariance  A(tj) 
(under  nominal,  hypothesized  conditions).  The  chi-square  random  variable  A!{tk)  is 
given  by 

^(ik)  =  (2-45) 

j=k-N+l 

with  N  being  the  size  of  a  sliding  detection  window.  Notice  that  the  system  dynamics 
are  not  included  in  Equation  (2.45)  and  that  only  one  failure  hypothesis  is  available. 
This  agrees  with  the  discussion  in  Section  1.5. 2. 3  which  stated  that  the  chi-square 
test  is  very  simple  and  can  only  be  used  to  detect,  not  identify,  failures.  A  detection 
rule  based  on  an  established  threshold  e  would  be 


N{tk)  >  e  Declare  FAILURE 

A{tk)  <  c  Declare  NO  FAILURE 


(2.46) 
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2.4  Summary 

This  chapter  provides  the  theoretical  basis  for  the  remaining  chapters.  MMAE 
will  be  used  for  FDI  in  this  work.  The  GLR  and  chi-square  algorithms  are  discussed 
here  because  they  will  be  used  as  a  benchmark  in  comparing  the  MMAE  FDI  re¬ 
sults  to  those  achieved  by  Vasquez  [41,42]  using  more  traditional  GLR/chi-square 
FDI.  The  theoretical  developments  of  this  chapter  closely  follow  those  presented  by 
Vasquez  [41,42]  and  Nielsen  [31]. 
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3.  Methodology 

3.1  Overview 

This  chapter  introduces  the  computer  software  and  simulation  techniques  used 
in  pursuing  this  research  effort.  An  overall  description  of  the  integrated  system  is 
given,  followed  by  detailed  state  and  measurement  models  for  each  of  the  navigation 
subsystems  used.  The  specific  multiple  model  filter  structure  used  to  detect  failures 
is  shown.  Finally,  the  simulation  approach  and  research  goals  are  reviewed  prior  to 
the  presentation  of  the  results  in  Chapter  4. 

3.2  Overall  System  Description 

The  main  elements  of  the  PLS  being  tested  for  failures  are  the  INS  and  the 
GPS  or  DGPS.  The  barometric  altimeter,  radar  altimeter,  and  ground-based  pseu- 
dolite  also  provide  measurements  to  the  Kalman  filter.  The  following  measurements 
are  available:  four  satellite  vehicle  (SV)  pseudoranges,  altitude  from  the  baromet¬ 
ric  altimeter,  one  surveyed-point  (pseudolite)  range  measurement  (optionally),  and 
height  above  ground  level  from  the  radar  altimeter  (optionally).  The  truth  model 
used  to  represent  the  real  world  consists  of  62  states,  while  the  Kalman  filter  model 
is  made  up  of  13  states. 

A  block  diagram  representing  the  system  PLS  configuration  is  shown  in  Fig¬ 
ure  3.1.  The  true  aircraft  position  is  generated  by  the  trajectory  profile  generator 
PROFGEN  [27]  and  provided  to  each  navigation  system.  The  GPS  satellite  vehicle 
(SV)  positions  are  given  by  actual  satellite  data  recorded  on  4  May  1991  and  are 
combined  with  the  true  aircraft  position  to  obtain  true  ranges,  which  are  modified 
with  noise  to  provide  pseudoranges  measurements  for  use  by  the  GPS.  Each  navi¬ 
gation  system  generates  perturbations  from  the  true  range  and  the  final  difference 
measurements  are  then  formed  by  subtracting  the  GPS  measured  ranges  from  their 
corresponding  INS- calculated  ranges.  The  EKF  equations  propagate  the  PLS  error 
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Feedforward  corrections 
(Indirect  method) 


Figure  3.1  Overall  PLS  Block  Diagram 

states  and  use  the  measurements  to  calibrate  its  state  estimates.  Finally,  these  state 
estimates  are  used  to  correct  the  INS-indicated  position  at  each  sample  time. 


3.3  PLS  Component  Model  Descriptions 

The  truth  model  consists  of  39  INS  states,  30  GPS  or  22  DGPS  states,  and  a 
single  pseudolite  state  (optionally).  The  filter  model  consists  of  11  INS  states  and 
two  GPS  or  DGPS  states.  The  following  sections  will  provide  the  details  of  these 
models  and  the  bases  for  their  selection. 

3.3.1  INS  Models. 

3. 3. 1.1  The  INS  Truth  Model.  This  section  presents  the  truth  model 
used  for  the  INS.  The  INS  is  a  strapped-down  wander  azimuth  system  based  on 
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the  Litton  LN-93.  The  manufacturer,  Litton  developed  a  93-state  error  model  [16] 
describing  the  error  characteristics  of  the  LN-93.  The  error  states  used  in  the  full 
model  may  be  separated  into  6  categories. 


SxfSxf SxJ SxJSxf Sxg 


where  is  a  93  x  1  column  vector  and 


(3.1) 


^Xi  represents  the  “general”  error  vector  containing  13  position,  velocity,  attitude, 
and  vertical  channel  errors;  the  first  nine  states  are  those  of  the  standard  Pinson 
model  of  INS  error  characteristics. 

Sx2  consists  of  16  gyro,  accelerometer,  and  baro-altimeter  exponentially  time-correlated 
errors,  and  “trend”  states.  These  states  are  modeled  as  first  order  Gauss- 
Markov  processes  in  the  truth  model. 

^X3  represents  gyro  bias  errors.  These  18  states  are  modeled  as  random  constants 
in  the  truth  model. 

5x4  is  composed  of  accelerometer  bias  error  states.  These  22  states  are  modeled  in 
the  same  manner  as  the  gyro  bias  states. 

5x5  depicts  accelerometer  and  initial  thermal  transients.  The  6  thermal  transient 
states  are  first-order  Gauss-Markov  processes  in  the  truth  model. 

5x6  models  gyro  compliance  errors.  These  18  error  states  are  modeled  as  biases  in 
the  truth  model. 
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The  truth  model  state  space  differential  equation  is  given  by 


A  full  description  of  the  submatrices  for  this  error  model  differential  equation  is  given 
in  Appendix  B.  This  93-state  error  model  is  a  highly  accurate  LN-93  representation, 
but  the  high  dimensionality  of  the  state  equation  makes  the  model  prohibitively 
CPU-intensive  (computationally,  and  in  terms  of  storage)  for  projects  examining  a 
large  number  of  problem  variations.  The  work  of  Negast  at  AFIT  addressed  the 
reduction  of  the  INS  error-state  model  [30]  while  preserving  enough  fidelity  to  be 
considered  a  viable  truth  model. 

The  reduced-order  model  to  be  used  as  the  truth  model  in  this  research  is 
defined  in  Equation  (3.3): 


This  model  was  shown  by  Negast  [30]  to  be  sufficient  to  represent  the  93-state  model 
accurately.  Note  that  the  submatrix  indices  used  in  representing  the  39-state  model 
are  not  identical  to  those  used  in  outlining  the  93-state  INS  error  model.  The 
relationship  between  the  two  models  is  given  in  Tables  A. 5  and  A. 6  of  Appendix  A. 
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3.3. 1.2  The  INS  Filter  Design  Model.  The  INS  filter  design  model  is 
the  model  which  would  be  used  by  the  EKF  operating  on  an  aircraft  using  the  PLS 
described  here.  The  limited  computational  power  available  and  the  requirement  for 
real-time  processing  motivates  making  the  dimension  of  the  filter  model  as  small  as 
possible.  The  INS  filter  model  is  comprised  of  11  states  (the  first  nine  being  the 
standard  Pinson  error  model  states):  6  misalignment  errors,  3  velocity  biases,  and 
2  states  for  barometric  stabilization.  Table  A. 9  shows  the  11  INS  states  used  in  the 
filter  model. 

3. 3. 1.3  The  INS  Measurement  Model.  The  only  measurement  model 
associated  with  INS  is  that  for  barometric  altimeter  aiding.  The  altimeter  aiding  is 
used  to  compensate  for  the  instability  inherent  in  the  vertical  channel  of  the  INS. 
The  altimeter  output  AltBaro  is  modeled  as  the  sum  of  the  true  altitude  ht,  the  total 
error  in  the  barometric  altimeter  Shs,  and  a  random  measurement  noise  v.  Similarly, 
the  INS  calculated  altitude  Altjj<[s  is  the  sum  of  the  true  altitude  and  the  INS  error 
in  vehicle  altitude  above  the  reference  ellipsoid,  8h.  A  difference  measurement  is 
used  to  eliminate  the  unknown  true  altitude,  ht,  resulting  in  Equation  (3.4): 

Sz  —  AltllSfS  —  Altsaro 

=  [ht  +  6h]-[ht  +  ShB-v]  (3.4) 

=  6h  —  Shs  A  V 

INS  error  in  vehicle  altitude  above  the  reference  ellipsoid,  8h,  and  total  barometric 
altimeter  correlated  error,  8hB,  are  states  10  and  11  in  the  11-  and  39-state  INS 
models. 

3.3.2  The  Radar  Altimeter  Model.  The  measurement  equation  of  the  radar 
altimeter  is  based  on  the  difference  between  the  INS-predicted  altitude  Altj^s  and 
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the  radar  altimeter  predicted  altitude  AltRait- 


8z  =  AltjtfS  —  Alt  Rah 

=  [h,  +  6h]-lh,-v]  (3.5) 

=  Sh  +  V 

The  errors  in  the  radar  altimeter  are  modeled  as  white  noise  with  no  time-correlated 
component.  This  may  be  a  rather  crude  model,  but  should  be  sufficient  to  demon¬ 
strate  performance  trends.  Note  that  no  additional  states  are  required  with  the 
addition  of  this  radar  altimeter  model. 

The  radar  altimeter  measurement  noise  covariance  RRait  is  a  function  of  aircraft 
altitude  above  ground  level  (AGL)  and  will  be  the  same  in  the  truth  and  filter  models. 
The  radar  altimeter  noise  covariance  from  [11]  is  altitude-dependent  and  is  given  by 

RRait  =  {[0.01]^  *[AGUrue?}  + 0.25  ft^  (3.6) 

3.3.3  GPS  Models.  The  GPS  (also  DGPS)  generates  user  position  based 
on  “known”  ranges  to  satellites  at  “known”  positions.  The  satellites  themselves 
transmit  their  position  in  space  (in  the  form  of  ephemeris  data)  as  accurately  as 
it  is  known  and  the  exact  time  (also  a  best  estimate)  at  which  the  transmission  is 
sent.  The  actual  range  information  is  calculated  based  on  knowledge  of  the  satellite 
position  and  the  finite  propagation  speed  of  the  electromagnetic  radiation  emitted 
from  the  satellite. 

3.3.3. 1  The  GPS  Truth  Model.  The  GPS  model  used  in  this  work 
was  developed  by  past  researchers  at  AFIT  [10,30,41,42].  The  dynamics  and  mea¬ 
surement  equations  for  the  full  30-state  truth  model  are  presented  in  this  section 
and  a  tabular  listing  of  the  states  is  shown  in  Table  A. 7  of  Appendix  A. 
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Five  types  of  error  sources  are  modeled  in  the  GPS  truth  state  equations.  Two 
of  the  five  error  sources  are  insignificant  when  differential  corrections  are  applied 
(DGPS).  The  first  error  type,  user  clock  error,  is  common  to  all  SV’s.  The  remaining 
four  error  types  are  unique  to  each  SV.  The  first  two  states  represent  user  clock  errors 
and  are  modeled  as; 


■ 

‘ 

' 

XUclkt 

0  1 

XUclkt 

>  = 

< 

XUclkar 

0  0 

^  XUclkir 

(3.7) 


where 

^Ucikt,  —  range  equivalent  of  user  clock  bias 
xuclkdr  —  velocity  equivalent  of  user  clock  drift 
The  initial  state  estimates  and  covariances  for  these  states  were  chosen  to  be  consis¬ 
tent  with  previous  AFIT  research  [3,10,30,41,42]  and  are: 

XUclk,,{to) 

< 

XUclkdri^o) 

and 

^Udkt,Uclkar{'^o)  = 

Because  these  error  sources  are  a  function  of  the  user  equipment,  they  are  common 
to  all  the  SV’s.  Recall  that  each  of  the  remaining  error  types  is  specific  to  each  SV, 
denoted  by  a  subscript  j. 

The  second  error  type  is  the  code  loop  error  dRdj.  The  code  loop  is  part  of  the 
user  equipment  shared  by  all  the  SV’s,  but  its  error  magnitude  is  relative  to  each 
SV.  The  work  of  Negast  [30]  shows  that  this  error  source  may  be  disregarded  in  the 
DGPS  model.  The  third  GPS  error  type  is  the  result  of  atmospheric  interference 
with  the  EM  signals  broadcast  by  each  SV,  specifically,  ionospheric  and  tropospheric 
delay,  SRionj,  and  SRtropj-  The  code  loop  error,  tropospheric  delay,  and  ionospheric 


9.0  X  10^^  ft^  0 

0  9.0  X  10^°  /tVsec" 


(3.9) 


^  — 

\ - 

o 

o 

1 _ 

(3.8) 
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delay  are  all  modeled  as  first-order  Markov  processes  with  time  constants  shown 
in  Equation  (3.10).  All  three  are  driven  by  zero-mean  white  Gaussian  noise  with 
strength  shown  in  Equation  (3.13).  The  fourth  error  source  is  due  to  inaccuracies 
of  the  clocks  on  board  the  individual  SV’s,  SRsdkj-  Like  code  loop  error,  this  error 
source  is  also  removed  when  differential  corrections  are  applied.  The  final  GPS  error 
source  is  based  on  line-of-sight  errors  between  the  SV’s  and  the  receiver,  8xsj ,  Sysj , 
and  8zsj.  The  DGPS  models  for  these  states  are  shown  later  in  Equations  (3.24)  - 
(3.27).  Note  that,  if  DGPS  were  used  exclusively  in  this  research,  the  states  for  code 
loop  error  and  satellite  clock  error  could  be  removed  completely.  The  modifications 
required  for  DGPS  are  summarized  in  Section  3. 3. 3.4. 


8Rcij 

-1 

0 

0 

8Rtropj 

0 

1 

500 

0 

8Rionj 

0 

0 

1 

1500 

^Rsclkj 

>  = 

0 

0 

0 

8Xsj 

0 

0 

0 

0 

0 

0 

8zsj 

\  J  J 

0 

0 

0 

with  initial  covariances  given  by 


0  0  0  0 

8Rdj 

'  •< 

Wd 

0  0  0  0 

8  Rtf  op  j 

Wtfop 

0  0  0  0 

8Rionj 

Wion 

0  0  0  0 

8Rsdkj 

>  +  < 

0 

0  0  0  0 

8x,- 

0 

0  0  0  0 

8ys, 

0 

0  0  0  0 

0 

(3.10) 


Pops  = 


0.25  ft^  0 

0  1.0  ft^ 
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0  0 
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0 
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(3.11) 
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and  noise  means  and  covariances  given  by 


£'[wGP5(i)]  =  0  (3.12) 

0.5  0  0  0  0  0  0 

0  0.004  0  0  0  0  0 

0  0  0.004  0  0  0  0 

E[wGPs{t)wGps{t  +  T)]  =  0  0  00000  ft'^/sec-S{T)  (3.13) 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

0  0  0  0  0  0  0 

The  full  30-state  GPS  dynamics  matrix  is  not  shown  explicitly  but  may  be  easily 
constructed  by  augmenting  Equation  (3.7)  and  four  copies  (one  for  each  SV)  of 
Equation  (3.10). 

3. 3. 3. 2  The  GPS  Filter  Design  Model.  Research  has  shown  [26,30] 
that  the  two  user  clock  error  states  provide  a  sufficient  filter  model  for  GPS.  The 
primary  argument  is  that  the  errors  modeled  for  the  28  other  GPS  (20  other  DGPS) 
states  (assuming  four  SV’s)  are  small  when  compared  to  the  user  clock  errors  which 
are  common  to  all  SV’s.  By  increasing  the  dynamics  driving  noise  and  re-tuning  the 
filter,  the  overall  performance  of  the  integrated  navigation  system  can  be  maintained. 
The  GPS  filter  model  is  given  by  Equation  (3.7)  plus  noise: 


iuclh 

XUclkdr 


Oil  XUclkt 

0  0  I  Xuclka^ 


+ 


'^clkdr 


(3.14) 


Because  simulations  were  only  performed  using  DGPS,  and  not  GPS,  no  experimentally- 
determined  Q  tuning  values  for  GPS  are  shown. 


S.3.3.3  The  GPS  Measurement  Model.  The  pseudorange  measure¬ 
ments  available  to  the  GPS  receiver  are  the  sum  of  the  true  range,  several  error 
sources,  and  a  random  noise; 

RopSj  —  Rtj  T  hRcioopj  T  hR^j-opj  T  hRionj  T  hRscikj  T  hRxjcik  (3.15) 

or,  after  differential  corrections  are  applied: 

RopSj  —  Rtj  ~i~  h R-ij.Qp^  h R^Qjij  -j-  V j  (3.16) 


where 


RcpSj  —  GPS  pseudorange  measurement,  from  SVj  to  user 

Rtj  =  true  range,  from  SVj  to  user 

hRcioopj  =  range  error  due  to  code  loop  error 

^Rtropj  =  range  error  due  to  tropospheric  delay 

6Rionj  =  range  error  due  to  ionospheric  delay 

SRscikj  =  range  error  due  to  SVj  clock  error 

hRudk  =  range  error  due  to  user  clock  error 

Vj  =  zero  —  mean  white  Gaussian  measurement  noise 


Because  Rt  is  not  available  to  the  filter,  a  substitution  will  be  make  to  eliminate  this 
term.  First,  the  satellite  position  vector  X5  and  the  user  position  vector  X?/  are 
defined  as: 


f  ' 

e 

Xg 

Xu  =  < 

Vu 

>  ,  Xs  =  < 

Vs 

Zu 

\  . 

Zs 

(3.17) 


where  the  superscript  e  denotes  coordinates  in  the  earth-centered  earth-fixed  (ECEF) 
frame.  The  pseudorange  from  the  user  to  the  satellites  calculated  by  the  INS,  Rins, 
is  the  difference  between  the  PLS-calculated  user  position,  Xy,  and  the  satellite 
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position  given  by  the  ephemeris  data,  X^: 


'  \ 

e 

e 

Xu 

Xs 

Rins  —  |Xf7  —  Xs\  — 

< 

yu 

.  —  . 

ys 

Zu 

Zs 

^  .. 

An  equivalent  form  of  Equation  (3.18)  is: 

Rins  =  -  ^sY  +  {yu  -  ysY  +  {su  -  zs)^ 


(3.18) 


(3.19) 


With  perturbations  representing  errors  in  Xjy  and  Xs,  Equation  (3.19)  can  be  written 
in  terms  of  the  true  range  and  a  truncated  first-order  Taylor  series: 


Rins  =  Rt^ 


9Xs 

,  9fljM's(Xg,Xu-) 
9Xi7 


(Xs,Xt;)„ 


SXs 


(Xs,Xc;)  nom 


SXu 


(3.20) 


The  solution  for  Rins  is  found  by  substituting  Equation  (3.19)  into  Equation  (3.20) 
and  evaluating  the  partial  derivatives  to  get: 


Rins 


(3.21) 


Finally,  the  truth  model  GPS  pseudorange  difference  measurement  is  given  as: 


8z 


Rins  —  Rgps 


+ 


•  8zu 

•  8zs 


8Rcioop  8R-iTop  dRi/Qii^^  -|-  v 


(3.22) 


The  user  position  errors  in  Equation  (3.22)  can  be  derived  from  the  first  three  (posi¬ 
tion  error)  states  of  the  filter  or  truth  model  using  an  orthogonal  transformation  [2]. 
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The  filter  design  measurement  model  for  the  GPS  measurement  does  not  con¬ 
tain  terms  for  the  errors  due  to  code  loop  variations,  atmospheric  delays,  satellite 
clock  deviations,  or  errors  in  ephemeris-given  satellite  position.  The  filter  GPS  mea¬ 
surement  model  can  be  written  as: 


6z 


Xs  -  Xu 


•  6xu  - 


ys  -  yu 

_  1-R/iVsl  . 


•  Syu  - 


Zs  -  Zu 


•  6zu  —  SRudk  +  V  (3.23) 


3. 3. 3. 4  The  DGPS  Truth  Model.  Differential  GPS  (DGPS)  is  mod¬ 
elled  very  similarly  to  GPS.  The  justifications  for  the  differences  were  given  above 
in  Sections  3.3.3. 1  and  3. 3. 3. 3  at  the  points  during  the  GPS  development  where  the 
differences  were  relevant.  This  section  will  present  the  mathematical  specifics  of 
DGPS  as  they  differ  from  those  of  GPS.  A  tabular  listing  of  the  22  DGPS  states  is 
shown  in  Table  A. 8  of  Appendix  A. 

The  error  sources  for  DGPS  are  identical  to  those  for  GPS,  with  the  noted 
exceptions  that  the  code  loop  and  satellite  clock  errors  may  be  disregarded  when 
using  differential  corrections.  The  differential  equation  for  the  DGPS  error  states 
becomes: 


0 

1 

1500 

0 

0 

0 


0  0  0  SRtropj 

0  0  0  hRionj 

0  0  0  ^  6xsj 

0  0  0 

0  0  0  6zs- 


(3.24) 
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with  initial  covariances  (note  the  differences  from  Equation  (3.11))  given  by 


1.0  ft^  Q  0  0  0 

0  1.0  0  0  0 

Pdgps  =  0  0  .35  0  0  (3.25) 

0  0  0  .35  ft‘^  0 

0  0  0  0  .35  ft^ 

and  noise  means  and  strengths  (note  the  differences  between  Equation  (3.13))  given 

by 

-E'[wx)GPs(f)]  =  0  (3.26) 

0.001  0  0  0  0 

0  0.0004  0  0  0 

E[wDGPs{t)^DGPs{t  +  t)]  =  0  0  0  0  0  ft'^lsec- 6{t)  (3.27) 

0  0  0  0  0 

0  0  0  0  0 

3.3.5.5  The  DGPS  Filter  Design  Model.  The  filter  design  model  for 
DGPS  is  identical  to  that  for  GPS.  The  only  difference  lies  in  the  gains  used  to 
tune  the  filter  when  using  DGPS.  The  filter  Q  values  used  to  tune  the  filters  in  this 
research  may  be  found  in  Appendix  C. 

3. 3. 3. 6  The  DGPS  Measurement  Model.  After  differential  corrections 
are  applied,  the  measurement  equation  is 

PdGPSj  —  Etj  T  ^Titropj  T  hRionj  T  hRiJclk  Pj  (3.28) 

as  noted  in  Section  3. 3. 3. 3.  After  noting  the  removal  of  the  code  loop  and  satellite 
clock  error  sources  in  Equation  (3.15),  leaving  Equation  (3.28),  the  remainder  of  the 
DGPS  measurement  model  development  is  identical  to  that  for  GPS. 


3-13 


3.3.4  The  Pseudolite  Model.  Pseudolite  measurements  are  treated  as  mea¬ 
surements  from  a  5*^  GPS  satellite,  with  the  following  three  exceptions  in  the  truth 
model. 

•  Transmissions  from  a  pseudolite  do  not  pass  through  the  ionosphere,  so  there 
is  no  ionospheric  delay  error  term  for  a  pseudolite  measurement.  The  tropo¬ 
spheric  delay  error  term  is  still  included. 

•  The  pseudolite  is  assumed  to  be  located  at  a  surveyed  position,  so  there  is  no 
uncertainty  in  the  “SV”  position  for  a  pseudolite  measurement. 

•  There  is  assumed  to  be  no  bias  or  drift  errors  in  the  “satellite”  clock  for  a 
pseudolite  measurement. 

As  may  be  inferred  from  the  list  above,  as  far  as  the  filter  is  concerned,  there  is  no 
difference  between  a  satellite  measurement  and  the  pseudolite  measurement. 

3.4  Failure  Models 

This  section  discusses  the  methods  used  to  model  failures  in  the  MMSOFE 
simulations.  For  each  failure  type  to  be  simulated,  the  corresponding  MMAE-based 
FDI  methods  used  to  detect  that  failures  are  also  discussed. 

3.4-1  Simulation  Failure  Models.  Seven  variations  of  three  different  failure 
types  are  modeled  in  simulation.  The  seven  failure  cases  are  presented  with  actual 
values  in  Table  3.1.  (The  contents  of  Table  3.1  will  be  discussed  in  detail  later  in 
this  section.)  A  short  description  of  each  failure  type  follows: 

1.  Interference/ Jamming  -  Interference  is  modeled  as  a  sudden  increase  in  the 
measurement  noise  associated  with  all  four  SV’s,  resulting  in  lower  carrier- 
to-noise  ratios,  CjNo,  in  the  GPS  receiver.  This  failure  is  induced  in  all  SV 
measurements  because  interference/jamming  is  assumed  to  occur  at  the  re¬ 
ceiver,  which  will  affect  all  four  (five,  if  the  pseudolite  is  included)  channels 
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simultaneously.  The  interference  noise,  is  added  to  the  truth  model  mea¬ 
surements  only.  Specific  magnitudes  of  interference  will  not  be  considered  as  in 
previous  work  at  AFIT  [41,42];  rather,  real-world  interference  will  be  allowed 
to  take  on  selected  values  within  the  interference  parameter  space  spanned  by 
the  MMAE  filter  bank.  Emphasis  will  be  placed  on  determining  (demonstrat¬ 
ing)  the  capability  of  MMAE  to  detect  and  identify  interference  failures  of 
unspecified  magnitude  quickly. 

GPS  jamming  is  used  to  refer  to  the  total  loss  of  useful  GPS  transmissions 
due  to  very  large  rf  interference.  A  GPS  jamming  failure  is  well- modelled 
(and  much  more  easily  modelled)  via  very  large  measurement  noise.  When 
the  MMAE  algorithm  detects  very  large  real-world  measurement  noise  R  then 
the  corresponding  measurements  will  be  very  lightly  weighted  by  the  elemental 
Kalman  filters;  the  effect  is  essentially  the  same  as  if  those  measurements  were 
never  received,  hence  the  use  of  the  term  “interference /jamming”  throughout 
this  report. 

2.  Spoofing  -  Spoofing  is  modeled  as  a  bias  added  to  the  measurements  associated 
with  all  GPS  SV’s.  The  addition  of  the  bias  (if  it  is  undetected)  has  the  effect  of 
placing  a  bias  on  the  position  solution  of  the  GPS  system.  Specific  magnitudes 
of  spoofing  failures  will  not  be  calculated  to  produce  some  effective  aircraft 
position  error,  according  to  a  single  strategy  of  spoofing,  as  the  emphasis  of 
this  work  is  an  assessment  of  the  EDI  capabilities  of  MMAE.  Instead,  a  range  of 
possible  spoof  magnitudes  will  be  investigated  to  exercise  the  MMAE  algorithm 
fully.  Two  models  of  spoofing  bias  addition  are  used  and  will  be  discussed 
shortly. 

In  this  research,  spoofing  is  modelled  as  a  uniform  bias  addition  to  each  of 
GPS  (and  pseudolite)  pseudorange.  Spoofing  failures  modelled  in  this  way 
result  in  (mainly)  vertical  changes  in  the  user’s  GPS  navigation  solution  (see 
Figure  3.2).  A  better  model  (more  representative  of  probable  malignant  spoof- 
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ing  techniques)  of  spoofing  would  be  the  addition  of  calculated  biases,  specific 
to  each  5y,  designed  to  yield  a  desired  net  change  in  the  user’s  navigation 
solution;  such  a  real-world  spoof  would  probably  be  designed  to  produce  a 
horizontal  shift  in  the  navigation  solution  so  that  unaffected  sensors  (baro  and 
radar  altimeters)  could  not  be  used  to  eliminate  the  error  over  several  measure¬ 
ment  updates.  Although  not  an  entirely  accurate  model  of  probable  real-world 
spoofing  changes,  the  uniform  bias  model  is  much  more  easily  implemented  in 
software  and  does  correctly  yields  cohesive  changes  in  the  user’s  GPS  naviga¬ 
tion  solution  which  effectively  test  the  spoofing  FDI  performance  of  the  MMAE 
algorithm. 

Line  of  constant  range  to  SV2  Line  of  constant  range  to  SVl 


Figure  3.2  Uniform-bias  spoofing  model:  vertical  effect 

(a)  Spoofing  Model  A  -  A  bias  is  added  suddenly  in  time  to  the  measure¬ 
ments  associated  with  each  SV,  with  the  net  effect  being  a  jump  in  the 
GPS-derived  aircraft  position.  More  sophisticated  spoofers  may  have  a 
fairly  accurate  estimation  of  the  GPS  measurements  being  received  by 
the  aircraft  and  can  discretely  add  a  small,  “undetectable”  step  bias  to 
these  measurements.  A  less  subtle  spooler  would  have  to  use  a  larger  bias 
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to  ensure  effectiveness  in  corrupting  the  measurements  while  running  the 
risk  of  being  detected. 

(b)  Spoofing  Model  B  -  A  more  intelligent  (and  considerably  more  difficult 
to  detect)  spoofing  failure  is  modeled  as  a  steadily  increasing  (ramp)  bias 
value  added  to  each  SV  pseudo-measurement.  Smaller  ramp  rates  are 
more  difficult  to  detect  than  large  rates  because  a  slow  change  may  be 
hidden  in  the  noise  expected  by  the  elemental  filters.  A  range  of  ramp 
rates  will  be  used  with  the  purpose  of  determining  the  detection  and  iso¬ 
lation  capabilities  of  MMAE.  Ramp  spoofing  will  be  termed  “intelligent” 
for  the  remainder  of  this  report. 

Table  3.1  summarizes  the  different  failures  that  are  simulated  in  this  research.  These 
failures  are  induced  on  each  of  the  PLS  hardware  configurations  under  test. 

Table  3.1  Failure  Types  and  Models  [41,42] 


L 

Fail  Type 

Description 

Fail  Method 

IQI 

No  Failure 

N/A 

R  =  Ro 

1 

Interference 

Jump  in 

Measurement  Noise 

Increase  R  from 

Ro  to  Ro  X  Magn^„j 

Step 

Spoofing 

Step  Bias  on  each 
Pseudorange 

Add  bias=Magn^py 
to  each  GPS  measurement 

Ramp 

Spoofing 

Ramp  Bias  on  each 
Pseudorange 

Add  bias=Magn^py  x  (t  —  to) 
to  each  GPS  measurement 

S.4-2  MMAE  Failure  Models.  A  bank  of  three  elemental  filters  in  an 
MMAE  structure  is  used  to  detect  and  isolate  interference  failures.  One  filter  is 
tuned  for  operation  when  no  interference  failure  has  occurred  (i.e.  measurement 
noise  R  =  Rq).  The  two  remaining  elemental  filters  are  each  tuned  for  best  operation 
with  increased  levels  of  measurement  noise  R.  An  increase  in  measurement  noise 
by  a  factor  of  2000  (i.e.  measurement  noise  R  =  Rox2000)  is  assumed  to  be  the 
highest  level  of  real-world  interference  that  will  be  encountered.  The  third  elemental 
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filter  models  this  highest  assumed  level  of  interference;  it  is  tuned  for  best  operation 
when  R  ==  Rox2000.  It  is  hypothesized  that  the  interference  elemental  filters  will 
operate  most  effectively  when  the  levels  of  their  respective  interference  assumptions 
are  separated  by  approximately  an  order  of  magnitude.  This  hypothesis  will  be 
verified  experimentally.  Based  on  this  hypothesis,  the  second  elemental  filter  is 
tuned  for  R  =  Rox200. 

MMAE-based  FDI  is  accomplished  via  the  information  contained  in  the  residu¬ 
als  of  the  elemental  filters.  The  residuals  produced  by  a  filter  with  an  accurate  model 
of  the  real  world  are  zero-mean,  white,  Gaussian,  and  have  covariance  HPH^  +  R. 
For  this  discussion,  we  will  term  these  residual  characteristics  “good” .  If  the  residu¬ 
als  of  any  elemental  filter  have  good  characteristics,  then  that  filter’s  hypothesis  of 
the  real  world  is  assumed  to  be  correct.  For  example,  if  the  R  =  Rox200  interference 
filter,  filter  two,  displays  good  residuals,  then  the  probability  of  model  correctness 
flows  to  filter  two  and  the  MMAF  algorithm  detects  a  real  world  interference  level 
200  times  Rq.  If  the  level  of  jamming  seen  by  the  filter  lies  between,  for  example, 
the  200xRo  and  2000xRo  levels  of  interference  modeled  by  elemental  filters  two  and 
three,  then  the  residuals  of  those  two  filters  will  both  display  some  good  character¬ 
istics.  In  this  case  the  probability  of  having  the  correct  model  will  be  shared  by 
models  two  and  three,  and  the  level  of  real  world  interference  is  isolated  by  blend¬ 
ing  the  measurement  noise  hypotheses  of  these  two  filters  based  on  the  computed 
probability  of  each  being  correct. 

Five  elemental  filters  are  used  to  detect  and  isolate  spoofing  failures  of  both 
step  and  ramp  types.  As  with  interference  failures,  the  first  elemental  filter  is  tuned 
for  best  operation  with  no  failure  induced.  The  second  elemental  filter  is  tuned  for 
best  operation  given  a  small  measurement  bias  added  to  the  SV  pseudoranges.  A 
third  elemental  filter  is  tuned  based  on  a  much  larger  measurement  bias  added  to 
the  SV  pseudoranges.  The  final  two  elemental  filters  model  correspondingly  small 
and  large  measurement  biases  of  a  negative  magnitude.  The  choice  of  this  spoofing 
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filter  bank  configuration  is  made  based  on  the  same  order-of-magnitude  separation 
arguments  used  to  hypothesize  the  configuration  of  the  interference  MMAE  filter 
bank.  In  each  case,  the  effectiveness  of  the  proposed  filter  bank  configurations  will 
be  evaluated  experimentally  and  changed  as  required  for  best  FDI  performance. 

For  step-type  spoofing  failures,  MMAF-based  FDI  is  conceptually  done  in  ex¬ 
actly  the  same  way  jamming  failures  are  detected  and  isolated  (although  the  next 
chapter  will  show  the  need  for,  and  implementation  of,  a  moving-bank  “pseudo¬ 
residual”  MMAF  to  handle  spoofs,  rather  than  a  standard  MMAF).  Probability 
flows  to  the  failure  models  (or  model,  in  the  case  of  an  exact  hypothesis  match)  that 
most  closely  match  the  real  world  spoofing  condition.  The  level  of  real  world  spoof¬ 
ing  is  estimated  by  blending  the  spoofing  hypotheses  of  the  elemental  filters  based 
on  the  computed  probability  associated  with  each  filter.  Because  of  the  blending 
capability  of  MMAF,  there  is  no  need  to  hypothesize  ramp  bias  failures.  As  the 
value  of  the  ramp  increases,  the  hypothesis  probability  should  gradually  flow  from 
one  model  to  the  next.  The  resulting  blended  estimate  of  the  magnitude  of  the  real 
world  spoofing  failure  will  follow  the  growth  of  the  actual  ramp  failure. 

The  final  steps  of  this  research  will  involve  identification  of  jamming  or  spoof¬ 
ing  failures  using  all  of  the  eight  elemental  filters  described  above  (one  nominal, 
three  jamming  strengths,  and  four  spoofing  biases).  These  simulations  will  establish 
the  performance  of  MMAF  FDI  under  the  operational-like  assumption  that  either 
jamming  or  spoofing  is  possible.  The  greatest  difficulty  for  the  FDI  algorithm  un¬ 
der  this  assumption  will  be  disambiguating  between  the  failure  types.  Performance 
analyses  will  argue  for  a  two-tiered  algorithm,  one  seeking  spoofs  and  the  second 
seeking  interference/jamming,  rather  than  a  single  MMAE  attempting  to  detect  all 
failure  modes  simultaneously. 


3-19 


3.5  Simulation  Software 

Multimode  Simulation  for  Optimal  Filter  Evaluation  (MSOFE)  is  a  general- 
purpose,  multimode  simulation  program  for  designing  integrated  systems  that  em¬ 
ploy  optimal  (Kalman)  filtering  techniques  and  for  evaluating  their  performance  [29]. 
The  general-purpose  construction  of  MSOFE  allows  its  application  to  a  wide  variety 
of  user-specific  problems  with  a  minimal  amount  of  new  software  development.  The 
United  States  Air  Force  uses  MSOFE  for  the  validation  of  all  systems  that  use  op¬ 
timal  filtering  techniques.  MSOFE  provides  Monte  Carlo  and  covariance  simulation 
modes. 

Multiple  Model  Simulation  for  Optimal  Filter  Evaluation  (MMSOFE)  was  de¬ 
veloped  by  Nielsen  [31,32]  at  the  Air  Force  Institute  of  Technology  to  support  the 
analysis  of  systems  using  a  multiple  model  adaptive  filter  structure.  MMSOFE  is 
written  as  an  extension  to  MSOFE  and  is  based  on  the  same  core  code.  The  MM¬ 
SOFE  program  propagates  multiple  filters  forward  in  parallel  while  performing  the 
hypothesis  probability  and  blending  calculations  required  for  MMAE  and  other  mul¬ 
tiple  model  algorithms.  The  Monte  Carlo  simulation  mode  of  MMSOFE  is  used  in 
all  phases  of  the  work. 

3.6  Summary 

This  chapter  has  presented  the  models  used  for  simulating  the  PLS,  includ¬ 
ing  the  navigation  component  models  and  the  failure  implementation  and  detection 
models.  The  failure  detection  and  simulation  methods  used  in  this  research  have  also 
been  discussed.  Results  and  analysis  of  these  simulations  are  presented  in  Chapter  4. 
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4-  Results  and  Analysis 


4.1  Introduction 

This  section  reviews  the  failure  models  and  navigation  component  configura¬ 
tions  used  in  this  research,  along  with  a  brief  discussion  of  why  each  is  of  interest.  A 
standard  MMAE  algorithm  will  be  used  subsequently  to  address  detection  of  jam¬ 
ming/interference  (increased  measurement  noise),  whereas  a  moving-bank  MMAE 
will  be  shown  to  be  more  appropriate  to  detect  and  compensate  for  the  onset  of 
spoofing  (measurement  bias). 

4.1.1  Failure  Models.  Jamming  is  the  total  loss  of  GPS  satellite  trans¬ 
missions  due  to  heavy  rf  interference.  With  the  loss  of  GPS  signals,  navigation  is 
totally  dependent  on  the  remaining  sensors,  and  for  this  study,  almost  entirely  based 
on  the  onboard  INS.  No  INS  currently  in  production,  or  even  realistically  foreseen, 
has  error  characteristics  small  enough  to  be  the  sole  primary  sensor  in  a  PLS.  While 
it  would  be  of  interest  to  observe  the  inflight  (not  landing)  navigation  performance 
of  an  MMAE  filter  bank  specifically  modelling  this  failure,  that  extension  is  beyond 
the  scope  of  this  work.  The  problem  of  jamming  was  briefly  considered  to  confirm 
the  above  presumptions.  These  results  are  presented  in  Section  4.2.2. 

RF  interference  not  strong  enough  to  cause  total  loss  of  the  GPS  transmissions 
is  modelled  as  increased  measurement  noise.  Interference  models  low-power  jamming 
from  hostile  sources,  as  well  as  navigation  in  areas  of  abnormally  high  rf  activity,  as 
in  the  immediate  area  of  television  or  radio  microwave  relay  stations.  Additionally, 
as  antenna  and  signal  processing  technology  improves,  devices  which  now  result  in  a 
total  loss  of  GPS  transmissions  may,  in  the  future,  only  have  the  effect  of  increasing 
the  measurement  noise  associated  with  the  GPS  receiver.  It  is  anticipated  that  the 
MMAE  algorithm  will  do  a  good  job  of  detection  and  isolation  against  this  failure 
type. 
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Spoofing  of  a  GPS  receiver  occurs  when  a  hostile  source  presents  GPS-like 
transmissions  to  the  receiver  which  are  slightly  stronger  than  the  real  GPS  signals 
and  which  will  give  an  incorrect  navigation  solution.  Intelligent  spoofers  would  be 
able  to  present  GPS-like  signals  giving  a  navigation  solution  with  very  little  initial 
offset  from  the  real  GPS  navigation  solution.  Once  the  receiver  accepts  these  false 
signals  they  are  “walked  off”  the  real  solution.  This  is  modelled  as  a  ramp  offset 
from  the  true  GPS  solution.  Less  intelligent  spoofers  would  present  signals  with  a 
very  large,  easily  identifiable  step  offset.  The  onset  of  unintelligent  spoofing  will  be 
seen  as  a  large  spike  in  the  measurement  residuals  within  the  elemental  filters  of  the 
MMAE  (or  within  a  single  non-adaptive  Kalman  filter  if  one  were  used  rather  than 
an  MMAE  algorithm).  Identification  is  easily  done  based  on  that  spike.  Intelligent 
spoofing  is  potentially  much  more  difficult  to  detect.  For  ramp  offsets  which  grow 
slowly  enough  to  be  hidden  in  the  noise,  there  are  no  means  available  for  detection. 

Realistically  speaking,  intelligent  spoofing  would  be  tremendously  difficult  to 
implement.  Prior  to,  and  over  the  duration  of,  the  spoof,  the  spoofer  would  have  to 
know  the  position  of  the  aircraft  with  very  nearly  the  same  precision  as  the  aircraft’s 
onboard  filters.  This  is  required  so  that  the  spoofer  could  maintain  a  slowly  growing 
offset  in  the  false  GPS  solution.  It  is  expected  that  the  MMAE  algorithm  will 
effectively  detect  and  isolate  spoofing  with  a  step  onset  or  a  ramp  onset  with  a  slope 
large  enough  to  be  distinguished  from  the  noise. 

4.1.2  Navigation  Component  Configurations.  The  four  navigation  com¬ 
ponent  configurations  used  are  shown  in  Table  4.1.  Case  one  represents  the  best 
configuration  available.  Cases  two  and  three  show,  respectively,  the  effect  on  EDI 
performance  of  eliminating  the  pseudolite  and  radar  altimeter  sensors.  Case  four 
shows  the  effect  of  a  much  lower  quality  INS,  which  might  be  representative  of  INS’s 
available  to  civil  aviation  groups. 
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Table  4.1  Navigation  Component  Cases 


Case 

INS  Type 

GPS  Type 

Altimeter  Aiding 

1 

0.4nm/hr 

DGPS  and 
Pseudolite 

Baro  and 
Radar  Alt 

2 

0.4nm/hr 

DGPS 

Baro  and 
Radar  Alt 

3 

0.4nm/hr 

DGPS 

Baro 

4 

4.0nm/hr 

DGPS 

Baro 

4-2  FDI  and  Navigation  Results 

4.2.1  Summary  Plots.  Two  types  of  plot  summaries  will  be  used  to  display 
the  simulation  results  of  this  study.  Each  simulation  was  performed  over  the  interval 
3810  to  3910  seconds,  which  represents  the  final  100  seconds  before  touchdown  in 
the  flight  profile  used  in  this  study.  A  lower  bound  probability  value  of  0.01  was  used 
in  all  cases  shown.  The  value  of  0.01  was  experimentally  determined;  smaller  lower 
bounds  slowed  the  response  of  the  ph's  to  onsets  of  their  corresponding  failures  (see 
Section  2.3.1).  Larger  lower  bound  values  did  not  provide  any  improvement  in  the 
Pk  response  time  and  reduced  the  total  probability  available  for  assignment  to  the 
filter  with  the  currently  correct  failure  hypothesis. 

Figure  4.1  is  representative  of  the  first  type  of  summary  plot.  The  first  subplot 
shows  the  magnitude  time  history  of  the  induced  failure,  indicated  by  dark  x’s.  The 
second  trace  on  the  first  subplot  shows  the  MMAE  filter  bank’s  blended  estimate  of 
the  real-world  failure  magnitude.  The  y-axis  of  this  first  subplot  is  the  magnitude  of 
the  total  variance  multiplier  (times  the  no-failure  noise  variance)  present/detected 
at  any  time  during  the  simulation.  Recall  that  interference  is  modelled  as  a  large 
increase  in  the  measurement  noise  variance.  The  remaining  subplots  show  the  15- 
run  mean  of  the  probability  (pk)  time  histories  for  each  of  the  elemental  filters.  The 
15-run  standard  deviation  of  each  of  the  pfc’s  is  also  shown.  The  y-axis  labels  of 
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Blended  Jamming  Solution  &  Elemental  Filter  Probabilities 


Figure  4.1  FDI  Performance  and  Elemental  Filter  Probabilities,  Interference  Fail¬ 
ures,  Three  Elemental  Filters  -  Nav  Case  1 
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these  probability  plots  describe  the  failure  assumptions  of  that  elemental  filter.  The 
reader  should  note  that,  as  is  expected,  the  blended  estimate  of  the  failure  magnitude 
visually  corresponds  to  the  sum  of  the  plotted  pk  values  times  their  respective  failure 
assumptions  as  indicated  in  the  y-axis  labels. 

Figure  4.2  is  representative  of  the  second  type  of  summary  plot.  Three  subplots 
will  always  be  included,  showing  the  latitude,  longitude,  and  altitude  errors  (in  feet) 
of  the  MMAE  blended  navigation  solution.  The  solid  trace  represents  the  blended 
filter-predicted  sigma  bound  for  that  variable.  The  dotted  line  shows  the  actual 
15-run  mean  of  the  same  variable  and  the  dashed  line  shows  the  15-run  mean  ±  one 
standard  deviation. 

4.2.2  Jamming.  Total  GPS  failure  is  modelled  in  this  work  by  zeroing  the 
H  (measurement)  matrix  entries  corresponding  to  the  GPS  measurements.  Better 
performance  might  be  realized  if  the  MMAE  bank  were  to  include  a  filter  only  mod¬ 
elling  INS  and  altimeter  states,  as  opposed  to  the  current  model  including  the  GPS 
states  and  using  the  zero  H  matrix  entries.  This  extension  was  not  a  primary  focus 
of  this  work  because  DGPS  jamming  is  well  modelled  (and  more  easily  simulated) 
as  very  strong  interference.  When  a  filter  assumes  a  very  large  noise  variance  R 
associated  with  some  measurements,  then  those  measurements  are  essentially  dis¬ 
carded  by  the  filter  when  it  calculates  its  state  updates.  Telling  the  filter  (via  large 
R  values)  not  to  use  some  measurements  has  the  same  effect  as  not  including  those 
measurements  at  all. 

4.2.3  Interference.  MMAE  EDI  was,  as  expected,  effective  against  interfer¬ 
ence  failures.  However,  as  can  be  seen  in  Figure  4.1,  the  filter  bank  has  some  difficulty 
in  blending  its  two  noise  hypotheses.  With  the  aid  of  hindsight,  it  could  have  been 
expected  that  the  MMAE  algorithm  would  be  prone  to  bouncing  from  hypothesis  to 
hypothesis  without  much  blending,  and  to  have  a  tendency  to  choose  the  model  with 
the  larger  assumed  noise  level.  As  the  real  world  noise  level  becomes  even  slightly 
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feet 


Figure  4.2  Navigation  Performance,  Interference  Failures  -  Nav  Case  1 
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greater  than  that  assumed  by  one  filter,  that  filter’s  residuals  very  quickly  look  bad 
(i.e.,  its  [r^(tj)A^^(ti)rfc(tj)]  becomes  a  large  value)  because  the  measurements  are 
consistently  violating  its  assumptions.  Although  the  measurement  noise  variance 
is  much  less  than  it  expects,  the  filter  assuming  a  larger  noise  value  sees  measure¬ 
ments  that  do  consistently  fall  within  its  expected  variance  (its  [r^(fi)A^^(ti)rfc(t,)] 
is  considerably  smaller).  Probability  quickly  flows  to  the  filter  with  the  larger  noise 
variance  assumption.  The  state  estimation  provided  by  the  three-filter  (assuming  lx, 
200x,  and  2000x  the  original  tuned  R  values)  interference  bank  is  shown  in  Figure  4.2. 
The  2000x  variance  level  was  chosen  as  the  upper-bound  of  possible  real-world  in¬ 
terference.  It  was  hypothesized  that  the  elemental  filters  should  be  separated  by 
approximately  an  order  of  magnitude  to  prevent  ambiguity.  The  200x  variance  filter 
was  selected  on  this  basis.  Filters  with  still  smaller  noise  variances  hypotheses  (20x 
for  example)  were  not  included  because  noise  variances  smaller  than  about  lOOx  have 
virtually  no  effect  on  the  navigation  system’s  performance.  A  larger  upper  bound 
on  the  assumed  real-world  noise  variance  might  just  have  easily  been  used,  resulting 
in  a  different  set  of  elemental  filters.  It  is  supposed  that  the  FDI  performance  of  the 
MMAE  algorithm,  as  shown  in  the  simulation  results,  would  not  differ  dramatically 
with  such  a  change.  As  can  be  seen  in  Figures  4.1  and  4.2,  the  estimation  perfor¬ 
mance  of  this  filter  bank  is  quite  good,  and  very  effectively  increases  the  filter  bank’s 
predicted  bounds  on  the  state  estimation  error  standard  deviations  and  so  prevents 
the  unacceptable  performance  or  filter  divergence  that  might  otherwise  result. 

Figures  4.3  and  4.4  show  the  performance  of  a  single  Kalman  filter  subjected  to 
interference  and  spoofing  failures,  respectively.  The  unacceptable  performance  seen 
in  these  figures  (note  the  scale  changes  from  Figure  4.2)  provides  both  a  benchmark 
for  performance  and  a  strong  motivation  for  adaptive  filtering  via  MMAE  techniques. 

Figure  4.5  shows  an  attempt  to  achieve  better  estimation  blending  by  including 
two  more  elemental  filters  with  interference  assumptions  in  the  same  range  spanned 
by  the  original  three  filters,  with  the  elemental  filters  now  assuming  lx,  lOOx,  300x, 
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Figure  4.3  Navigation  Performance,  Interference  Failures,  single  Kalman  filter  - 
Nav  Case  1 
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1600x,  and  2000x  the  original  tuned  R  values.  The  values  of  these  noise  variance 
assumptions  were  chosen  somewhat  arbitrarily,  with  the  purpose  of  testing  the  ear¬ 
lier  hypothesis  that  an  order  of  magnitude  should  separate  the  variance  assumptions 
of  the  elemental  filters.  As  might  be  expected  (again  with  the  aid  of  hindsight)  the 
filters  close  together  at  the  upper  end  of  the  interference  range  tend  to  confuse  each 
other  as  their  assumed  noise  variances  differ  by  only  twenty  percent.  Additional 
blending  can  be  seen  at  low  interference  values,  where  the  elemental  filter’s  R  values 
are  still  separated  by  substantial  percentages,  however,  as  is  seen  in  the  figure,  that 
this  blending  takes  place  does  not  guarantee  that  better  estimation  performance  is 
realized  in  the  low  true  R  range.  The  probability  plot  summary  of  this  bank  is 
included  only  to  justify  the  selection  of  the  three-filter  interference  bank  for  final 
implementation  and  so  the  state  estimation  of  the  five-filter  interference  bank  is  not 
shown.  These  results  suggest  the  use  of  an  MMAE  filter  bank  spanning  the  interfer¬ 
ence  range  of  interest  and  composed  of  elemental  filters  separated  by  approximately 
an  order  of  magnitude  in  their  parameter  assumptions. 

4.2.4  Spoofing.  A  five-filter  (zero  measurement  bias,  ±15  foot  bias,  and 
±240  foot  bias)  MMAE  bank  was  used  for  measurement  bias  detection.  Assumed 
bias  values  of  ±15  feet  were  included  to  allow  the  bank  to  fine  tune  its  bias  estimation. 
Filters  separated  by  15  feet  are  close  enough  in  parameter  space  to  provide  clear 
blending  and  are  far  enough  away  to  avoid  ambiguity  between  the  filter  models. 
The  bias  assumption  of  the  ±240  foot  bias  filter  pair  was  determined  as  follows.  It 
was  found  experimentally  that  an  elemental  filter  was  able  to  detect  actual  biases 
within  150  feet  of  its  own  assumed  bias,  while  biases  further  removed  than  150  feet 
resulted  in  zero  (lower  bound  probability)  pk  values  for  those  elemental  filters.  It 
was  decided  that,  for  good  blending,  one  half  of  this  detection  range  should  overlap 
the  detection  range  of  the  next  closest  filter,  e.g.,  the  225  foot  bias  assumption 
difference  between  the  two  positive  filters  at  15  feet  and  240  feet  is  obtained  as 
150  feet(range)  ±  150  feet(range)  -  75  feet(overlap).  These  filter  assignments  are. 
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of  course,  problem-specific.  The  reader  will  note  that  each  of  the  elemental  filters 
hypothesizes  a  constant  bias  offset.  Ramp  spoof  offsets  are  estimated  over  time  by 
blending  the  constant  bias  hypotheses  of  the  same  elemental  filters. 

The  following  paragraph  discusses  the  summary  plots  of  two  spoofing  simula¬ 
tions  in  order  to  motivate  (1)  the  development  of  a  modified  MMAE  algorithm  to 
detect  measurement  offset  failures,  and  (2)  the  use  of  moving-bank  MMAE  in  con¬ 
junction  with  this  modification.  The  questions  raised  in  the  next  paragraph  will  be 
answered  in  Section  4.2.4. 1,  which  develops  the  theory  of  the  MMAE  modification 
mentioned  above.  Throughout  the  remainder  of  this  report,  this  modified  MMAE 
algorithm  will  be  referred  to  as  “pseudo-residual”  MMAE  (PRMMAE)  for  reasons 
discussed  in  Section  4.2.4. 1. 

Figure  4.6  shows  the  spoofing  detection  performance  of  the  presumedly  stan¬ 
dard  MMAE  algorithm  (it  was  later  realized  that  the  MMAE  algorithm  actually 
implemented  in  the  MMSOFE  software  was  a  non-moving-bank  pseudo-residual 
MMAE,  to  be  described  in  the  next  section).  In  order  to  observe  the  detection 
behavior  of  the  algorithm  clearly,  each  real-world  spoofing  level  change  was  chosen 
to  match  exactly  the  spoofing  level  modelled  by  one  of  the  elemental  filters.  The 
no-spoof  filter  correctly  acquires  the  GPS  satellites  and  initializes  the  state  estima¬ 
tion.  At  the  onset  of  the  first  spoofing  block,  the  probability  of  the  elemental  filter 
assuming  the  negative  of  the  real  world  spoof  value  spikes  for  one  sample  period 
(see  Section  4. 2. 4.1).  After  this  one  sample  period  of  useful  information  from  the 
MMAE  filter  bank,  the  probabilities  of  all  of  the  elemental  filters  become  equal  (see 
Section  4. 2. 4.1),  suggesting  that  some  kind  of  re-centering  of  the  filter  bank  at  the 
initial  time  when  information  is  available  may  prove  useful.  Good  state  estimation  is 
not  recovered  even  when  the  real  world  bias  is  removed.  While  searching  for  a  means 
of  improving  the  meager  detection  performance  shown  in  Figure  4.6,  it  was  discovered 
that  the  the  MMAE  algorithm  had  been  incorrectly  implemented  with  respect  to  the 
addition  of  measurement  biases  (“pseudo-residuals”  versus  true  residuals  were  used 
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Figure  4.6  Spoofing  FDI  Performance  and  Elemental  Filter  Probabilities,  Non- 
Moving  Bank  PRMMAE  -  Nav  Case  1 
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Figure  4.7  Spoofing  FDI  Performance  and  Elemental  Filter  Probabilities,  Standard 
MMAE  -  Nav  Case  1 
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to  form  [rjA'^^rk]  for  the  MMAE’s  probability  calculations;  see  Section  4. 2. 4.1). 
Figure  4.7  shows  the  detection  performance  of  the  corrected  (see  Section  4.2.4. 1) 
MMAE  algorithm.  This  performance  clearly  shows  that  correct  MMAE  is  unable  to 
detect  measurement  bias  offsets.  The  total  probability  is  equally  divided  among  the 
five  elemental  filters  because  all  of  the  filters  have  indistinguishable  residuals.  The 
“ELF1”-“ELF5”  labels  seen  in  Figure  4.7  were  used  early  in  this  research  to  denote 
ELemental  Filters  one  through  five,  and  correspond  to  bias  offset  assumptions  of  0, 
+350,  -350,  +700,  and  -700  feet,  respectively.  State  estimation  performance  plots 
are  not  included  for  Figures  4.6  and  4.7  because  these  cases  are  included  only  to 
motivate  the  development  of  the  pseudo-residual  MMAE  (PRMMAE)  algorithm. 

4. 2. 4-1  Pseudo-Residual  MMAE  Theory.  The  incorrect  MMAE  al¬ 
gorithm  produces  good,  potentially  exploitable  information  for  one  sample  period 
while  the  correct  MMAE  algorithm  does  not.  The  mathematical  reasons  for  this 
unexpected  result  follow. 

At  all  times 

^true  ~  Hfrue^true  +  ^true  +  '^true  (4.1) 

Let  us  assume  two  filter  models,  (1)  assuming  no  measurement  bias,  i.e.,  that 
hirue  =  0,  and  (2)  assuming  a  positive  measurement  bias,  i.e.,  that  htme  —  bi,  where 
bi  will,  at  least  initially,  be  the  bias  actually  simulated  in  the  real  world.  These  two 
filters  will  have,  as  their  best  prediction  z  of  the  measurement  before  it  arrives 


(1)  Zi  =  Hxi 

(2)  Z2  ==  Hxj  +  bi 


(4.2) 


which  will  produce  the  following  update  equations  used  to  generate  x+  at  each  sample 
period: 

(1)  Xf  =  X7+Ki[Zt,„e-Zi] 

(2)  xj"  =  X2  +K2[ztr,ie  -Z2] 
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(1)  xf  =  Xi  +Ki[Zfr„e  -  Hxi] 

(4.4) 

(2)  Sit  =  X2  +K2[Zir„e-(Hx2  +bi)] 

After  allowing  both  filters  to  run  to  steady  state,  each  filter  will  modify  its  state 
estimates  so  that  the  residual  is  zero-mean,  i.e.  E[ztrue  —  z]  =  0.  In  fact,  it  was 
observed  that  the  filter’s  primarily  altered  their  user  clock  bias  state  to  yield  such 
zero- mean  residuals.  This  has  the  following  effect: 


(1)  E[HX]^  —  ^truet^irue  —  O) 


—  H<,.ueXf7. 


(2)  E[Hx2  ]  —  '^truei^true  —  0)  Ej  —  'H.true^true  Ei 

At  the  onset  of  a  spoof,  htrue  becomes  non-zero,  htme  —  Ei  for  example.  At  the  next 
measurement  update,  the  true  residuals  are  (discounting  noise) 

(1)  [l^true  Zi]  =  [(HjrueXtrue  T  Ej)  HXj  ]  r  a  r>\ 

(4.6) 

(2)  [Zirue  Z2]  =  [(HtrueXtrue  "f"  Ej)  (HX2  -b  Ej)] 

and  the  pseudo-residuals,  namely  [ztme  —  Hx^]  for  all  k,  are  (discounting  noise) 

(1)  [Zirue  HXj  ]  —  [(HtrueXtrue  T  Ei)  HXj  ] 

(2)  [Zirue  HX2  ]  —  \{^H.true^true  T  El)  HX2  ] 

The  true  residuals,  using  E[Hx^]  from  Equation  (4.5),  become 


(1)  [H-lrueXt  rue  "f*  bi  —  bj 

(2)  [H^T’^igXtT’ue  ”1“  bi  ((Hfrue^true  bi)  -|-  bi)]  —  bi 

and  the  residuals,  using  E[Hx^]  from  Equation  (4.5),  become 

(1)  rue  T  El  HirweXtrue]  —  Ei 

(2)  [HjrueXiT-ue  T  El  (HtrueXjj-ue  Ei)]  =  2Ei 
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In  writing  Equations  (4.1)  to  (4.9),  it  is  assumed  that  true  residuals  are  used 
to  update  the  elemental  filters;  what  is  subject  to  consideration  here  is  whether  the 
true  residuals  or  the  pseudo- residuals  are  more  useful  in  forming  for  the 

MMAE’s  probability  computations.  As  may  be  seen  from  Equation  (4.8),  the  true 
residuals  from  each  elemental  filter  do  not  reveal  the  real  world  measurement  bias 
and  are  indistinguishable  from  one  another  (see  Figure  4.6).  It  may  be  deduced 
from  Equation  (4.9)  that  the  filter  assuming  the  negative  of  the  actual  bias  will 
show  nearly  zero-mean  pseudo-residuals  at  the  measurement  update  immediately 
following  the  spoof  onset  (recall  Figure  4.6).  Hence,  good  identification  may  be 
achieved  by  using  the  pseudo-residuals  to  form  In  the  above  discussion, 

had  a  filter  assumed  a  bias  of  — bi,  then  in  steady  state  its  Hx"  would  have  become 
HtrueXirue  +  bj.  The  pseudo-residuals  of  such  a  filter  would  be  roughly  zero- mean 
at  the  measurement  update  following  the  spoof  onset.  The  information  provided  by 
this  zero-mean  measurement  pseudo-residual  is  used  to  isolate  the  actual  bias.  The 
true  residuals  must  be  used  (as  in  a  single  Kalman  filter  or  regular  MMAE)  to  update 
the  elemental  filters. 

4.2. 4-2  Pseudo-Residual  MMAE  Performance.  Examination  of  Fig¬ 
ure  4.6  in  the  light  of  the  previous  section  suggests  that  good,  reliable  information 
does  exist  with  which  to  identify  the  onset  of  a  spoof  (isolation  is  also  possible  if  an 
individual  filter  happens  to  assume  exact  value  of  bias  offset  that  is  applied),  but  it 
is  only  visible  for  a  single  measurement  update.  The  moving-bank  MMAE  algorithm 
introduced  here  is  useful  in  exploiting  this  information.  The  reader  should  note  that, 
up  until  this  point,  we  have  limited  our  discussion  to  a  fixed-bank  MMAE. 

To  get  around  the  quirk  of  only  having  the  decision  information  available  for 
one  update  cycle,  one  only  has  to  accept  the  single  sample  period  of  information 
as  an  identification  of  the  spoof,  subsequently  estimate  the  spoof,  then  go  back 
and  reprocess  the  last  measurement  assuming  that  the  estimated  spoof  is  present  in 
the  real  world  [the  entire  filter  bank  is  moved  to  the  neighborhood  of  the  estimated 
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spoof).  Filter  bank  movement  is  accomplished  by  subtracting  hesUmated  from  the  true 
measurements  before  they  enter  the  MMAE  algorithm,  rather  than  adding  hestimated 
to  each  filter’s  Zjt,  since  it  is  computationally  more  efficient.  This  process  is  repeated 
for  that  single  update  time  until  the  new  spoof  bias  value  is  completely  identified  and 
the  bank  is  recentered.  After  this  correction,  each  of  the  elemental  filters  in  the  bank 
steps  forward  into  the  next  propagation  cycle  without  knowing  it  ever  experienced  a 
spoof.  This  process  can  be  observed  in  the  pk  plots  of  simulations  involving  spoofing 
elemental  filters.  In  the  first,  third,  and  fifth  subplots  of  Figure  4.10  (to  be  seen 
and  described  in  detail  later  in  this  section,  associated  with  ramp  rather  than  step 
spoofs,  but  called  out  here  because  it  clearly  shows  the  desired  phenomenon),  for 
example,  it  can  be  seen  that  when  the  real-world  spoof  bias  changes  by  -1-15  feet/sec 
(the  first  spoofing  ramp),  the  pk  associated  with  the  -15  foot  bias  elemental  filter 
spikes  for  one  second.  The  bank  is  moved  based  on  this  spike,  and  in  the  next  second 
the  probability  returns  to  the  no-bias  elemental  filter.  Throughout  the  simulation 
depicted  in  Figure  4.10  (and  each  spoofing  simulation),  it  may  be  seen  that  the  pk 
value  corresponding  to  the  no-bias  elemental  filter  displays  such  downward  spikes 
each  time  a  change  in  the  real-world  spoof  bias  value  is  detected.  These  downward 
spikes  depict  the  information  given  by  the  first  measurement  update  which  caused 
the  MMAE  filter  bank  to  be  moved;  in  reality,  due  to  the  measurement  reprocessing 
described  above,  the  no-bias  elemental  filter’s  pk  value  never  moves  from  its  near¬ 
unity  value. 

Figure  4.8  shows  the  detection  and  isolation  performance  of  the  PRMMAE  al¬ 
gorithm  in  the  moving-bank  configuration.  The  reader  should  note  several  aspects  of 
Figure  4.8  which  are  consistent  with  the  description  of  the  moving-bank  PRMMAE 
algorithm’s  operation  as  described  so  far.  It  can  be  seen  that  the  probability  rests 
consistently  (excepting  the  downward  spikes  mentioned  above)  on  the  no-bias  ele¬ 
mental  filter  (fourth  sub-plot),  while  the  real-world  and  declared  spoof  values  range 
over  2000  feet  (first  sub-plot).  This  is  a  clear  indication  that  the  no-fail  filter  is  con- 
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sistently  recentered  on  the  real-world  spoof  value.  Each  time  the  real-world  spoof 
value  jumps  by  a  value  modelled  in  one  of  the  elemental  filters,  the  probability  spikes 
on  the  filter  assuming  the  negative  of  the  spoof  jump  value  (sub-plots  five-eight).  The 
third  sub-plot  (measurement  count)  shows  the  number  of  measurement  updates  that 
were  taken  to  identify  the  new  spoof  offset  value  fully  at  each  sample  period  before 
the  following  propagation  cycle  was  entered. 

Figure  4.8  and  the  corresponding  state  estimation  error  plots  of  Figure  4.9 
indicate  that  spoofing  jumps  are  identified  accurately  regardless  of  whether  their 
value  is  exactly  modelled  in  the  bank  of  elemental  filters.  The  PRMMAE  algorithm 
itself  is  unable  to  identify  spoofing  values  greatly  different  from  those  modelled  in 
its  bank.  The  difficulty  here  is  that  the  MMAE  residual  information  term,  r^A“^r, 
becomes  very  large  when  the  spoof  offset  is  numerically  greater  than  150  (feet) 
displaced  from  the  bias  assumption  of  a  modelled  filter.  Because  —1/2  times  this 
large  term  appears  in  an  exponential,  the  MMAE  conditional  hypothesis  calculation 
value  goes  to  zero  and  the  MMAE  algorithm  cannot  make  a  decision  about  which 
direction  to  move  the  bank  in  the  possible  failure  space.  The  first  proposed  solution 
to  this  problem  was  to  attempt  to  move  the  existing  bank  throughout  the  possible 
“spoof  offset”  parameter  space,  doing  measurement  updates  at  each  assumed  bias 
offset,  until  the  new  spoof  was  encountered.  As  it  turns  out,  this  is  not  necessary. 
Even  though  the  MMAE  calculations  using  the  term  become  useless,  the 

individual  measurement  residuals  associated  with  the  DGPS  measurements  in  each 
elemental  filter  allow  simple  isolation  of  the  spoof  magnitude.  These  residuals  are 
zero-mean  up  until  the  addition  of  the  spoof.  When  the  spoof  occurs,  it  shows  up 
directly  on  the  residuals.  Estimation  is  a  matter  of  simply  reading  the  number.  It 
can  be  seen  in  Figure  4.8  (sub-plots  1,  5-8)  that  large  changes  in  the  real-world  spoof 
value  correspond  to  probability  values  of  0.2  =  If  {^filters)  for  each  filter.  When 
all  filters  in  the  bank  have  equally  bad  residuals,  the  probability  is  equally  divided 
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Figure  4.9  Spoofing  (Steps)  Navigation  Performance,  Moving-Bank  PRMMAE  - 
Nav  Case  1 
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between  the  filters  and  the  individual  residual  terms  are  examined  for  the  spoof  offset 
as  described  above. 

Identification  of  spoofing  as  described  in  the  previous  paragraph  is  accom¬ 
plished  equally  well  with  only  two  symmetrical  bias-assuming  filters  displaced  about 
a  filter  assuming  no  bias  offset,  and  not  the  four  symmetrical  filters  that  were  used 
earlier  in  this  research.  The  three-filter  spoofing  bank  is  used  to  identify  the  spoof 
onset  and  to  fine-tune  the  spoof  estimation  once  the  first  guess  is  made.  Based  on 
the  experience  gained  during  this  research,  and  given  the  methods  finally  used  to 
detect  measurement  bias  failures,  it  may  be  possible  to  identify  the  spoof  with  only 
one  filter.  In  fact,  Vasquez’  [41,42]  GLR  methods  of  spoof  identification,  on  first 
inspection,  appear  to  be  mathematically  quite  similar  to  the  methods  used  in  this 
research  but  using  a  single  Kalman  filter.  Bank  movement  and  fine-tuning  the  spoof 
isolation  are  only  easily  performed  via  MMAE  blending  using  the  no-bias  and  two 
symmetrical-bias  filters. 

Figure  4.10  shows  the  FBI  performance  of  the  moving-bank  PRMMAE  algo¬ 
rithm  in  the  face  of  ramped  (intelligent)  spoofing  offsets.  The  first  ten  real-world 
spoofing  ramp  segments  (sub-plot  one)  have  slopes  of  30,  24,  20,  16,  12,  10,  8,  6,  4, 
and  2  feet  per  second,  respectively.  It  can  be  seen  that  no  significant  identification 
error  occurs  until  the  8  ft/sec  spoofing  segment  (note  the  dotted  meanicr  trajectories 
on  the  first  subplot).  The  smaller  spoofing  ramps  do  cause  considerable  confusion  to 
the  algorithm,  although  the  state  estimation  still  does  not  go  too  far  awry,  as  seen 
in  Figure  4.11.  It  can  be  concluded  that,  given  the  navigation  component  configu¬ 
ration  of  case  one,  spoofing  steps  and  spoofing  ramps  as  small  as  10  ft/sec  can  be 
identified  with  no  significant  error.  The  FBI  performance  of  cases  2-4  is  discussed  in 
Section  4.2.6.  The  blended  navigation  solution  provided  by  the  PRMMAE  algorithm 
is  quite  good  and  is  shown  in  Figures  4.9  and  4.11.  Comparison  of  these  results  with 
those  of  a  single  non-adaptive  Kalman  filter  (refer  to  Figure  4.4)  shows  the  dramatic 
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improvement  in  performance  given  by  the  moving-bank  PRMMAE  algorithm  used 
in  this  research. 

An  important  consideration  is  inherent  throughout  the  above  discussion.  The 
elemental  filters  of  the  MMAE  must  be  allowed  to  initialize  with  no  measurement 
bias  present  in  the  real  world.  What  is  actually  detected  and  tracked  is  the  offset  of 
the  measurement  bias  from  what  was  present  during  initialization.  This  assumption 
is  reasonable  because  GPS  initialization  will  presumably  be  done  before  take-off  from 
a  friendly  air  base.  Any  attempted  spoofing  there  would  be  detected  by  surveyed 
receivers  at  the  station,  and  shortly  removed. 

4.2.5  Detection  of  Both  Spoofing  and  Interference  in  a  Single  MMAE.  Fig¬ 
ure  4.12  shows  the  EDI  performance  of  a  filter  bank  containing  some  elemental  filters 
assuming  interference  and  some  assuming  spoofing  offsets.  Figure  4.13  shows  the  cor¬ 
responding  state  estimation  error.  The  five  modelled  filters  assume  (1)  zero  bias,  no 
interference  (lx  original  measurement  noise  variance),  (2)  zero  bias,  interference  rep¬ 
resented  by  200x  original  R,  (3)  zero  bias,  interference  modelled  by  2000x  original  R, 
(4)  -|-15  feet  bias,  no  interference,  and  (5)  -15  feet  bias,  no  interference.  The  inclu¬ 
sion  of  the  spoofing  filters  is  not  expected  to,  and  in  fact  does  not,  detract  from  the 
MMAE’s  performance  against  interference  because  increased  random  measurement 
noise  looks  nothing  like  the  uniform  measurement  bias  shown  by  a  spoof.  Detection 
of  interference  is  accomplished  with  no  degradation  of  performance  versus  the  case 
when  no  spoofing  filters  were  present.  However,  using  only  a  single  sample  period 
test,  there  is  no  possible  way  for  the  algorithm  to  distinguish  between  the  onset  of 
a  constant  measurement  bias  and  a  large  measurement  noise  value  superimposed  on 
the  unbiased  measurement.  In  an  MMAE  filter  bank  containing  both  interference 
and  spoof  elemental  filters,  the  spoof  elemental  filters  do  not  accept  probability,  at 
each  sample  time,  for  instances  of  random  interference  or  small  spoof  steps  because 
the  interference  elemental  filters  are  pre-disposed  to  declare  failures.  The  increased 
measurement  noise  R  assumption  of  the  interference  elemental  filters  causes  their 
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A  matrices  to  be  much  larger  than  the  A  matrices  associated  with  the  spoofing 
elemental  filters  (A  =  HPH^  +  R).  As  a  result,  the  elemental  filter  probability  cal¬ 
culation  exponential  {  — |r^A“^r}  is  a  negative  scalar  number  of  smaller  magnitude 
for  interference  filters,  because  of  their  larger  A  matrices,  than  for  spoofing  filters 
assuming  correct  bias  offsets,  even  when  the  spoofing  filters  have  somewhat  smaller 
residuals  r.  The  interference  elemental  filters  have  a  tendency  to  declare  failures 
in  this  case  because  they  are  (incorrectly)  more  heavily  weighted  than  the  spoofing 
elemental  filters. 

As  might  be  supposed,  identification  in  this  case  (spoof  and  interference  el¬ 
emental  filters)  is  possible  if  the  measurement  bias  added  is  of  greater  magnitude 
than  is  consistent  with  a  modelled  noise  variance  increase.  The  magnitudes  of  the 
first  nine  spoofing  steps  (sub-plot  one)  are  9000,  8000,  7000,  6000,  5000,  4000,  3000, 
2000,  and  1000  feet,  respectively  (offsets  are  bigger  for  this  test  because  smaller  spoof 
offsets,  as  in  the  earlier  simulations,  were  misidentified  as  noise).  Spoofing  identi¬ 
fication  is  done  correctly  until  the  2000  foot  spoofing  change.  Generally  speaking, 
interference  and  spoofing  cannot  be  distinguished  by  combining,  into  a  single  MMAE 
algorithm,  the  methods  used  in  this  research  to  identify  interference  and  spoofing 
separately. 

Future  researchers  may  find  great  success  identifying  both  interference  and 
spoofing  using  only  filters  with  spoofing  assumptions  and  extending  the  detection 
test  period  over  two  or  more  sample  periods.  If  the  residual  deviations  and  estimated 
spoof  value  display  a  constant  bias,  then  spoofing  may  be  identified.  If  residual 
deviations  and  estimated  spoof  value  instead  display  substantial  sample-to-sample 
changes,  then  increased  measurement  noise  variance  may  be  assumed  and  estimated 
(over  the  course  of  several  sample  periods  used  to  build  an  estimate  of  the  noise 
strength).  This  conjecture  of  spoofing  bank  behavior  when  subjected  to  interference 
failures  was  not  verified  in  this  research  due  to  a  shortcoming  of  the  simulation  tool 
MMSOFE.  MMSOFE  assigns  a  new  random  noise  realization  each  time  a  single 
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measurement  is  reprocessed,  making  it  impossible  for  the  spoof  filters  to  recenter  on 
the  original  noise-corrupted  measurement.  Note  that  this  noise  realization  concern 
does  imply  that  MMSOFE  was  also  incorrectly  assigning  a  new  measurement  noise 
realization  each  time  the  PRMMAE  algorithm  recentered  the  filter  bank  against  a 
real-world  spoof;  recentering  in  this  case  was  done  correctly  however  because  the 
magnitude  of  the  spoof  bias  b  dominated  relatively  small  magnitudes  of  the  incon¬ 
sistent  measurement  noise  realizations. 

One  can  envision  a  hierarchical  MMAE  structure  that  initially  looks  only  for 
spoofs  (using  a  moving-bank  PRMMAE  algorithm  structure),  but  if  the  estimated 
spoof  exhibits  large  sample-to-sample  changes,  the  elemental  filters  in  the  MMAE 
could  be  redefined  to  look  for  interference  instead  (via  a  non-moving-bank  standard 
MMAE  algorithm  structure).  After  the  noise  variance  is  estimated,  the  MMAE  el¬ 
emental  filters  can  be  returned  to  look  for  spoofs  only  (and  the  algorithm  returned 
to  a  moving-bank  PRMMAE  form),  but  now  with  each  elemental  filter  being  tuned 
for  the  correctly  estimated  measurement  noise  variance.  The  greatest  difficulty  in 
this  method  will  be  the  (possibly)  required  storage  of  the  state  estimation  during  the 
several  sample  periods  required  to  isolate  the  nature  of  the  failure.  Once  the  failure 
type  and  magnitude  is  identified,  the  MMAE  state  estimates  might  be  reprocessed 
over  that  period  of  time.  This  storage  and  reprocessing  may  or  may  not  be  necessary. 
That  will  need  to  be  determined  via  trial  and  error  to  see  whether  good  estimation 
is  maintained  with  the  simpler  method  of  performing  sample-to-sample  estimation 
and  not  reprocessing.  For  example,  it  may  be  acceptable  to  use  the  state  estimates 
produced  by  an  MMAE  looking  only  for  spoof  offsets  (with  its  moving-bank  struc¬ 
ture)  even  if  the  decision  is  made,  over  a  few  sample  periods,  that  interference  more 
accurately  models  the  current  real-world  failure.  Once  the  failure  is  isolated,  it  may 
be  required  to  reprocess  all  of  the  propagation  and  update  cycles  that  were  used  to 
isolate  the  failure,  using  the  newly  declared  bank  definition.  Alternatively,  a  parallel 
set  of  MMAE’s,  one  searching  for  spoofs  and  the  other  for  interference,  could  be 
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employed;  again  the  filter  seeking  only  spoofs  would  be  used  as  the  primary  algo¬ 
rithm,  but  if  the  decision  is  made  to  declare  an  interference  rather  than  a  spoofing 
failure,  the  alternate  MMAE  might  be  usable  without  delay  for  any  reprocessing  of 
the  propagation  and  update  cycles. 

The  use  of  a  jamming/spoofing  sequence  (described  shortly)  presents  a  very  dif¬ 
ficult  challenge  for  any  detection  algorithm.  A  hostile  enemy  employing  this  method 
would  heavily  jam  the  area,  causing  the  loss  of  the  GPS  satellite  signals.  The  en¬ 
emy  would  next  introduce  spoofed  GPS-like  signals,  then  remove  the  jamming  so 
that  the  GPS  receiver  re-acquires  the  spoofing  signals  rather  than  the  actual  GPS 
signals.  The  MMAE  EDI  algorithm  (or  any  other  algorithm)  does  not  maintain 
good  state  estimation  during  the  jamming  portion  of  this  failure.  When  the  jam¬ 
ming  is  removed,  the  filter’s  state  estimation  is  so  far  off  that  it  has  no  choice  but 
to  acquire  whatever  GPS  or  GPS-like  signal  is  present  in  the  real  world.  MMAE 
fails  at  detecting  the  spoofing  signal  following  the  jamming  phase.  The  best  that 
can  be  hoped  for  against  such  a  failure  is  that,  following  the  removal  of  jamming, 
spoofing  values  greater  than  the  possible  drift  of  the  onboard  INS  will  be  detected. 
Although  the  MMAE  EDI  would  have  difficulty  detecting  the  initial  offset  due  to 
spoof  once  the  jamming  were  removed,  it  could  compensate  for  step  or  ramp  spoof 
offsets  thereafter,  relative  to  that  initial  offset.  This  failure  sequence  was  not  tested 
during  this  research.  The  inclusion  of  an  elemental  filter  only  modelling  INS  and 
altimeter  states,  i.e.,  not  including  DGPS  states  rather  than  assuming  zero  H  matrix 
entries  corresponding  to  the  DGPS  measurements,  may  give  useful  information  in 
the  face  of  this  failure  sequence.  Using  this  filter,  at  re-acquisition,  any  spoofing 
signals  with  an  offset  larger  than  the  possible  INS  drift,  would  be  rejected  and  the 
spoof  identified. 

4.2.6  Performance  Comparisons:  Navigation  Cases  2-4.  This  section 
shows  the  failure  detection  and  isolation  performance  available  using  a  navigation 
component  suite  of  lower  quality  (and  less  cost)  than  that  used  earlier  in  this  chap- 
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ter  to  demonstrate  the  performance  of  the  MMAE  and  PRMMAE  algorithms  (see 
Table  4.1).  Unlike  the  plots  presented  previously  in  this  chapter,  which  are  grouped 
to  show  the  EDI  and  corresponding  navigation  performance  or  a  given  failure  case, 
the  plots  of  this  section  will  be  organized  to  emphasize  the  changes,  first,  of  EDI 
performance,  and  second,  of  navigation  performance,  with  the  range  of  navigation 
suites  given  in  Table  4.1.  Case  2  depicts  the  effect  of  removing  pseudolites,  case 
3  the  further  impact  of  removing  the  radar  altimeter,  and  case  4  the  compounded 
impact  of  using  a  lower-precision  INS.  The  EDI  performance  of  each  navigation  case 
against  interference  and  spoofing  step  failure  types  are  examined. 

Figures  4.14  through  4.17  show  the  interference  EDI  performance  of  navigation 
component  cases  one  through  four.  Note  that  the  speed  and  accuracy  of  detection 
and  isolation  of  interference  are  nearly  identical  for  each  of  the  four  cases. 

Figures  4.18  through  4.21  show  the  spoofing  EDI  performance  of  navigation 
component  cases  one  through  four.  It  can  be  seen  in  these  four  plots  that  spoofing 
detection  occurs  in  one  second  independent  of  the  navigation  suite  used.  The  degree 
of  misidentification  of  the  spoofing  magnitude  shown  in  the  second  subplot  of  Fig¬ 
ures  4.18  through  4.21  suggests  that  estimation  of  the  magnitude  of  spoofing  failures 
is  also  independent  of  the  navigation  components  used. 

Figures  4.22  through  4.25  show  the  navigation  performance  of  navigation  cases 
one  through  four  in  the  presence  of  interference  failures.  Examination  of  these  figures 
reveals  an  expected  degradation  of  the  navigation  performance  from  navigation  case 
one  to  two,  two  to  three,  and  three  to  four.  Comparison  of  Figure  4.23  (navigation 
case  two)  with  Figure  4.22  (navigation  case  one)  shows  the  effect  on  navigation 
performance  of  removing  the  ground-based  pseudolite  from  the  navigation  suite. 
Comparison  of  Figure  4.24  (navigation  case  three)  with  Figure  4.23  (navigation  case 
two)  shows  the  effect  on  navigation  performance  of  removing  the  radar  altimeter.  It 
can  be  seen  that  this  change  has  a  fairly  substantial  impact  on  the  precision  of  the 
aircraft  altitude  estimates,  as  anticipated.  Comparison  of  Figure  4.25  (navigation 
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Figure  4.14  FDI  Performance,  Interference  Failures  -  Nav  Case  1 
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Figure  4.15  FDI  Performance,  Interference  Failures  -  Nav  Case  2 
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Figure  4.16  FDI  Performance,  Interference  Failures  -  Nav  Case  3 
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Figure  4.17  FDI  Performance,  Interference  Failures  -  Nav  Case  4 
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Figure  4.18  FDI  Performance,  Spoofing  Step  Failures  -  Nav  Case  1 
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Figure  4.19  FDI  Performance,  Spoofing  Step  Failures  -  Nav  Case  2 


4-37 


•240  Bias  +240  Bias  -15  Bias  +15  Bias  No  Bias  MSMTCNT  Spoof  iD  Err 


3820  3 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 

10i - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 

H  C - [ - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 3 

I  r  K/ - — 1 1 - 1 1 - 1 1 — 1 1 - 1 1 - 1 1 - 1 1 - 1 1 - 1 1 — 1 1 1 1 - 1 1 - 1 1 — 1 1 - 1 1 - n 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 

- ^ ^ - ; - , - ^ ^ ^ ^ ^ ^ - 3 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 

- ^ ^ ^ ^ - , - , - ^ ^ ^ - n 


3810  3820  3830  3840  3850  3860  3870  3880  3890  3900  3910 

Time  (sec) 


Figure  4.20  FDI  Performance,  Spoofing  Step  Failures  -  Nav  Case  3 
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case  four)  with  Figure  4.24  (navigation  case  three)  shows  the  effect  on  navigation 
performance  of  using  a  navigation  suite  with  a  much  less  precise  INS.  Because  of  the 
consistently  good  FDI  performances  for  these  four  cases,  the  navigation  performance 
for  each  case  is  essentially  that  of  a  single  extended  Kalman  filter  correctly  designed 
and  tuned  for  the  navigation  system  components  and  artificially  informed  of  the 
actual  interference  noise  variance. 

Figures  4.26  through  4.29  show  the  navigation  performance  of  navigation  cases 
one  through  four  when  subjected  to  spoofing  step  failures.  Similar  observations  may 
be  made  about  the  relative  navigation  precision  offered  by  the  different  navigation 
suites. 


4.2.7  Comparison  to  GLR/Chi-Square  FDI.  The  FDI  performance  of  the 
MMAE  techniques  used  in  this  work  compare  very  favorably  with  the  GLR/chi- 
square  techniques  applied  to  a  similar  problem.  Due  to  time  constraints,  direct 
comparisons  of  MMAE  and  GLR/chi-square  using  the  detection  problems  of  this 
work  were  not  made;  however,  several  conclusions  may  be  drawn  based  on  the  similar 
FDI  study  performed  by  Vasquez  [41,42].  Vasquez’  study  was  similar  (dis-similar) 
in  the  following  ways:  Vasquez  used 

1.  essentially  the  same  INS  as  that  used  for  navigation  component  case  1  of  this 
study 

2.  GPS  rather  than  DGPS 

3.  no  radar  altimeter  and  no  pseudolites 

4.  a  six-transponder  range  range-rate  velocity  aiding  system 

The  GLR/chi-square  scheme  used  by  Vasquez  [41,42]  on  a  GPS-aided  inertial 
navigation  system  is  effective  at  detecting  interference  failures,  but  this  identification 
comes  after  a  (sometimes  large)  delay.  In  Vasquez’  [41,42]  work,  identification  of 
the  onset  of  interference  required  a  delay  of  2  seconds.  The  MMAE  algorithm  used 
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Figure  4.22  Navigation  Performance,  Interference  Failures  -  Nav  Case  1 
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Figure  4.23  Navigation  Performance,  Interference  Failures  -  Nav  Case  2 
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Figure  4.25  Navigation  Performance,  Interference  Failures  -  Nav  Case  4 
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Figure  4.26  Navigation  Performance,  Spoofing  Step  Failures  -  Nav  Case  1 
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Figure  4.27  Navigation  Performance,  Spoofing  Step  Failures  -  Nav  Case  2 
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Figure  4.28  Navigation  Performance,  Spoofing  Step  Failures  -  Nav  Case  3 
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Figure  4.29  Navigation  Performance,  Spoofing  Step  Failures  -  Nav  Case  4 
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in  this  work  identifies  and  estimates  the  magnitudes  of  interference  failures  with  a 
one  sample  period  delay  (one  second).  Identification  of  the  interference  magnitude 
is  not  attempted  using  the  GLR/chi-square  FDI  methods.  A  much  longer  delay, 
46  seconds  in  one  case,  is  required  for  the  GLR/chi-square  algorithm  to  return  to 
a  nominal  no-fail  declaration  from  a  large  interference.  MMAE’s  also  suffer  from  a 
longer  delay  in  returning  to  the  nominal  condition  from  a  large  failure,  compared 
to  the  time  required  to  move  either  from  a  nominal  condition  to  a  failure  condition 
or  from  one  failure  condition  to  another.  The  first  sub-plot  of  Figure  4.1  shows,  for 
example,  that  ten  seconds  are  required  for  the  MMAE  algorithm  implemented  in 
this  research  to  return  from  a  failed  to  a  no-fail  declaration.  Nevertheless,  this  is  far 
less  than  the  46  second  delay  of  the  GLR/chi-square  alternative. 

Spoofing  type  failures  (step  and  ramp)  were  detected/identified  less  effectively 
by  the  GLR/chi-square  detection  scheme  than  were  interference  failures.  Vasquez’  [41, 
42]  simulations  showed  that  spoofing  step  failures  with  magnitudes  as  low  as  50  feet 
were  detected  and  identified  (accurate  state  tracking  regained)  in  about  20  seconds. 
This  report  shows  that  the  PRMMAE  algorithm  detects  and  identifies  spoofing  steps 
as  small  as  15  feet  with  a  one  second  delay.  The  largest  (ramp)  spoofing  failure  sim¬ 
ulated  by  Vasquez  was  2  feet/sec  and  required  a  minimum  of  250  seconds  (chi-square 
test)  to  detect  this  failure.  The  smallest  spoofing  ramp  that  was  effectively  detected 
and  identified  by  PRMMAE  techniques  (in  one  second)  was  about  8  feet/sec;  how¬ 
ever,  it  can  seen  by  close  examination  of  Figure  4.10  that  a  spoofing  ramp  of  2 
feet  / sec  is  detected  within  about  5  seconds  even  though  reliable  isolation  at  this  low 
magnitude  is  prevented  by  system  noise. 

It  may  be  concluded  that  GLR/chi-square  failure  testing  has  a  significant  as¬ 
sociated  delay  when  compared  to  MMAE  failure  testing,  especially  in  the  face  of 
spoofing  (bias-like)  failures.  The  selection  of  one  or  the  other  of  these  FDI  schemes 
for  application  to  a  particular  problem  should  be  made  while  fully  taking  into  con¬ 
sideration  this  performance  difference. 
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4-3  Conclusions 

This  chapter  shows  the  development  of  moving-bank  pseudo-residual  MMAE 
identified  as  PRMMAE,  a  new  technique  for  the  identification  of  measurement  offset 
(bias  or  ramp  “spoof”)  failures.  PRMMAE  (for  spoofs)  and  standard  MMAE  (for 
interference)  are  used  to  detect  and  compensate  for  interference  and  spoofing  failures 
in  the  GPS  portion  of  a  navigation  configuration.  State  estimation  before,  during, 
and  after  these  failures  is  maintained.  These  results  show  that  a  precision  landing 
system  (PLS)  based  on  these  navigation  components  can  reliably  detect  degradation 
of  the  navigation  solution  due  to  external  rf  sources,  and  can  preserve  the  quality  of 
navigation  so  that  flight  and  some  categories  of  instrument  landings  may  safely  be 
continued. 
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5.  Conclusions  and  Recommendations 

Chapter  4  presented  and  analyzed  the  results  obtained  using  MMAE  and  the 
newly  developed  moving-bank  pseudo-residual  MMAE  (PRMMAE)  to  detect  and 
isolate  interference  and  spoofing  failures  in  a  DGPS-aided  inertial  system.  This 
chapter  briefly  summarizes  and  extends  the  conclusions  suggested  by  the  simulation 
results  and  analysis  discussed  at  length  in  Chapter  4. 

5.1  Introduction 

Much  recent  research  and  FAA  interest  has  been  directed  toward  the  develop¬ 
ment  of  a  system  to  aid  landing  navigation  to  replace  the  aging  instrument  land¬ 
ing  system  (ILS)  now  in  use.  It  is  widely  assumed  that  the  replacement  for  the 
ILS  will  be  based  on  the  global  positioning  system  (GPS).  Previous  research  at 
AFIT  has  resulted  in  the  development  of  a  DGPS-aided  (and  radar  altimeter-aided) 
INS-based  precision  landing  system  (PLS)  capable  of  meeting  the  FAA  precision 
requirements  for  instrument  landings.  The  susceptibility  of  DGPS  transmissions  to 
interference/jamming  and  spoofing  must  be  addressed  before  DGPS  may  be  safely 
used  as  a  major  component  of  such  a  safety-of-flight  critical  navigational  device. 

This  thesis  applies  multiple  model  adaptive  estimation  (MMAE)  techniques  to 
the  problem  of  detecting  and  identifying  interference/jamming  or  spoofing  failures 
in  the  DGPS  signal.  Interference/) amming  failures  are  modelled  as  significantly 
increased  measurement  noise  associated  with  the  DGPS  pseudorange  measurements. 
Spoofing  failures  are  modelled  as  constant-  (or  ramped-)  offset  values  added  directly 
to  the  DGPS  pseudorange  measurements. 

5.2  FDI  Conclusions 

The  poor  navigation  performance  of  a  non- adaptive  Kalman  filter  subjected  to 
interference/) amming  or  spoofing  motivated  the  investigation  of  MMAE  as  a  means 
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Figure  5.1  MMAE  and  Non-Adaptive  Kalman  Filter  Navigation  Performance  Com¬ 
parison,  Interference/ Jamming  Failures 

of  failure  detection  and  identification.  MMAE  is  very  effective  at  detecting  and 
isolating  interference/ jamming  failures.  The  navigation  performance  achieved  us¬ 
ing  MMAE  is  greatly  improved  from  the  non-adaptive  case.  The  left-hand  side  of 
Figure  5.1  shows  the  detection,  identification,  and  resulting  Mended  navigation  per¬ 
formance  of  the  MMAE  filter  bank  subjected  to  a  jamming  signal.  The  right-hand 
side  of  Figure  5.1  shows  the  performance  of  a  single  non-adaptive  Kalman  filter  with 
the  same  jamming  signal  applied.  (Note  the  scale  differences  of  the  two  altitude 
error  subplots).  Because  of  the  good  FDI  performance  of  the  MMAE  algorithm  (see 
Figure  5.1,  first  subplot,  left-hand  side),  the  navigation  performance  is  essentially 
that  of  a  single  extended  Kalman  filter  correctly  designed  and  tuned  for  the  navi¬ 
gation  system  components  and  artificially  informed  of  the  actual  interference  noise 
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variance.  Figure  5.1  pertains  to  the  case  of  a  high-precision  INS  aided  by  both  DGPS 
(with  pseudolite)  and  a  radar  altimeter,  but  the  conclusion  remains  valid  even  as 
pseudolite  and  radar  altimeter  measurements  are  removed,  and  the  high-precision 
INS  is  replaced  with  an  inexpensive  INS  with  accuracy  degraded  by  an  order  of 
magnitude. 

When  subjected  to  spoofing  (measurement  bias)  failures,  both  the  non-adaptive 
Kalman  filter  and  a  standard  MMAE  algorithm  (designed  to  seek  spoofing)  yield  ex¬ 
tremely  poor  navigation  solutions.  This  thesis  shows  the  development  of  a  moving- 
bank  pseudo-residual  MMAE  (PRMMAE)  algorithm  to  detect  and  identify  spoofing 
failures  in  the  DGPS  measurements.  The  PRMMAE  algorithm  uses  the  pseudo¬ 
residuals,  rather  than  the  conventional  residuals,  to  form  [r^A“^r]  for  use  in  the 
MMAE  hypothesis  probability  calculations.  Moving-bank  PRMMAE  is  very  effec¬ 
tive  at  detecting  and  identifying  spoof  bias  and  ramp  failures;  the  resulting  naviga¬ 
tion  performance  is  equivalent  to  that  of  a  single  extended  Kalman  filter  operating 
in  a  no-fail  environment.  Figure  5.2  shows  the  navigation  performance  of  the  PRM¬ 
MAE  algorithm  compared  to  that  of  a  non-adaptive  extended  Kalman  filter  when 
subjected  to  a  sequence  of  spoofing  bias  failures  (topmost  subplots). 

Detection  and  isolation  of  bias-like  (or,  even  worse,  ramp-like)  failures  on  the 
measurement  signal  is  a  difficult  task  for  any  FDI  algorithm,  including  MMAE.  Per¬ 
haps  the  most  significant  result  of  this  research  is  the  development  and  verification  of 
the  moving-bank  PRMMAE  algorithm  (see  Section  4.2.4. 1)  to  detect  and  isolate  such 
bias-like  or  ramp-like  failures.  It  is  hypothesized  that  moving-bank  PRMMAE  will 
prove  to  be  a  useful  and  generally  applicable  tool  for  the  detection  and  identification 
of  signal  bias-like  or  ramp-like  failures. 

Sections  4.2.3  and  4.2.4  show  the  FDI  effectiveness  of  standard  MMAE  and 
moving-bank  PRMMAE  against  DGPS  interference  and  spoofing  failures,  respec¬ 
tively.  It  may  be  generally  observed,  based  on  the  results  of  these  sections,  that 
detection  of  interference  or  spoofing  failures  (by  two  separate  algorithms  or,  presum- 
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PRMMAE  vs.  Spoofing  Single  Non-Adaptive  Filter  vs.  Spoofing 
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Figure  5.2  MMAE  and  Non-Adaptive  Kalman  Filter  Navigation  Performance  Com¬ 
parison,  Spoofing  Step  Failures 

ably,  by  two  algorithms  invoked  in  a  hierarchical  form,  but  not  by  a  single  filter  bank: 
see  Section  4.2.5)  is  accomplished  in  a  single  sample  period  (one  second).  Identifica¬ 
tion  of  correct  parameter  values  associated  with  these  failures,  in  general,  requires 
between  one  and  three  seconds.  This  performance  is  exceptional  when  compared  to 
a  GLR/chi-square  algorithm,  (arguably)  the  next  best  FDI  alternative  in  terms  of 
effectiveness  (see  Section  4.2.7). 
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5.S  Navigation  Conclusions 

Detection  and  identification  of  interference/ jamming  and  spoofing  failures  is  ac¬ 
complished  very  effectively  using  MMAE  and  moving-bank  PRMMAE,  as  discussed 
in  the  previous  section.  The  resulting  navigation  performance  of  the  integrated  sys¬ 
tem  with  these  failures  induced  is  only  slightly  degraded  from  that  attained  with  no 
failures,  rather  than  unacceptably  degraded  as  is  the  case  for  a  single  non-adaptive 
extended  Kalman  filter  as  conventionally  used  for  an  integration  filter  for  opera¬ 
tional  aided  inertial  navigation  systems  (compare  Figure  4.2  to  4.3  and  Figure  4.9 
to  4.4).  Table  5.1  shows  a  listing  of  the  navigation  component  cases  tested;  these 


Case 

INS  Type 

GPS  Type 

Altimeter  Aiding 

1 

0.4nm/hr 

DGPS  and 
Pseudolite 

Baro  and 
Radar  Alt 

2 

0.4nm/hr 

DGPS 

Baro  and 
Radar  Alt 

3 

0.4nm/hr 

DGPS 

Baro 

4 

4.0nm/hr 

DGPS 

Baro 

Table  5.1  Navigation  Component  Cases 


cases  were  selected  to  demonstrate  the  FDI  and  navigation  performance  impact  of 
eliminating  (or  degrading)  the  individual  navigation  sensors,  with  the  hope  of  ex¬ 
tending  the  results  of  this  research  to  potential  medium-  and  low-cost  commercial  or 
civilian  aviation  applications.  As  was  hoped,  the  failure-induced  degradation  of  nav¬ 
igation  performance  in  the  integrated  system,  based  on  the  MMAE  and  moving-bank 
PRMMAE,  for  all  navigation  component  cases  and  failures  tested,  is  not  substantial 
enough  to  prevent  sufficiently  accurate  navigation  to  achieve  mission  requirements 
in  other  (not  landing)  phases  of  flight. 

It  is  assumed  that  commercial  and  civilian  aircraft  will  be  much  more  likely 
to  face  interference/jamming  GPS  failures  (due  to  the  low  cost  of  implementing  an 
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interference  device)  than  spoofing  failures.  Of  particular  note,  then,  is  the  FDI  per¬ 
formance  of  each  navigation  component  case  against  interference /jamming  failures. 
Section  4.2.3  shows  that  the  interference/jamming  FDI  performance  of  navigation 
component  case  four  of  Table  5.1  is  essentially  identical  to  that  of  navigation  compo¬ 
nent  case  one.  Spoofing  identification  is  also  accomplished  very  effectively  by  all  four 
navigation  component  cases  (see  Section  4.2.4).  As  was  desired,  the  FDI  and  navi¬ 
gation  performance  of  these  cases  is  good  enough,  and  the  cost  of  these  lower-quality 
navigation  components  is  low  enough,  to  warrant  investigating  the  application  of  the 
PLS  to  commercial  and  civil  aviation  as  well  as  to  military  aircraft. 

5.4  Recommendations  for  Future  Research 

Future  researchers  may  find  good  results  by  pursuing: 

•  Hierarchical  structure  to  accomplish  FDI  of  jamming/interference  and  spoofing 
simultaneously. 

The  results  of  Chapter  4  show  that  a  composite  MMAE  filter  bank  (some  ele¬ 
mental  filters  looking  for  spoof-  and  some  for  interference-type  failures)  cannot 
disambiguate  between  interference/jamming  and  spoofing  failures.  One  can  en¬ 
vision  a  hierarchical  MMAE  structure  that  initially  looks  only  for  spoofs  (using 
a  moving-bank  PRMMAE  algorithm  structure),  but  if  the  estimated  spoof  ex¬ 
hibits  large  sample-to-sample  changes,  the  elemental  filters  in  the  MMAE  could 
be  redefined  to  look  for  interference  instead  (via  a  non-moving-bank  standard 
MMAE  algorithm  structure).  After  the  noise  variance  is  estimated,  the  MMAE 
elemental  filters  can  be  returned  to  look  for  spoofs  only  (and  the  algorithm  re¬ 
turned  to  a  moving-bank  PRMMAE  form),  but  now  with  each  elemental  filter 
being  tuned  for  the  correctly  estimated  measurement  noise  variance. 

•  Use  of  dither  signals  to  enhance  navigation  performance  of  the  integrated  sys¬ 
tem. 
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Much  research  has  demonstrated  the  usefulness  of  applying  a  known  dither 
signal  to  improve  the  identifiability  of  a  system  [8,19,38].  The  application  of 
such  a  dither  signal  has  the  effect  of  “shaking”  up  the  modes  of  the  system 
that  might  otherwise  be  difficult  to  observe  and  use  for  identification  purposes. 
While  the  FDI  and  navigation  performance  results  obtained  in  this  research 
are  quite  good,  no  attempt  was  made  to  apply  a  dither  signal  (or  even  to  corre¬ 
late  the  dynamics  of  the  flight  profile  used)  to  achieve  improved  performance. 
Adding  purposeful  input  dithers  to  enhance  failure  identifiability  may  warrant 
further  investigation,  particularly  in  phases  of  flight  prior  to  final  approach, 
flare,  and  landing. 
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Appendix  A.  Error  Model  State  Definitions 


This  appendix  contains  a  tabular  listing  of  the  93  Litton  INS,  39  reduced-order 
INS  truth  model,  30  GPS,  22  DGPS,  and  13  simulation  filter  error  states.  Note  that 
the  INS  model  is  shown  with  41  states.  This  is  because  the  two  barometric  aiding 
states  are  included  as  states  10  and  11  of  this  model.  The  71-state  (72-state  when  the 
pseudolite  is  included)  PLS  truth  model  states  are  defined  in  the  right-most  column 
of  Tables  A. 5  through  A. 8. 

A.l  Litton  LN-93  Error-States 

Tables  A.l  through  A. 4  list  the  LN-93  error  model  (93  states)  as  defined  in  the 
Litton  CDRL  [16].  Note  that  the  Litton  document  contains  several  errors  which  are 
corrected  in  these  tables  [30]. 

A. 2  Reduced  Order  INS  Truth  Model  States 

Tables  A. 5  and  A. 6  list  the  39  INS  states  (plus  two  barometric  aiding  states) 
used  in  all  the  full-order  models  in  this  thesis.  This  reduced  order  model  was  devel¬ 
oped  and  verified  by  Negast  [30].  States  12  and  13  of  the  PLS  are  defined  later  in 
Tables  A. 7  and  A. 8. 

A. 3  GPS  Error  States 

Table  A. 7  lists  the  GPS  error  states.  A  total  of  30  states  are  included  to  model 
the  error  characteristics  of  4  space  vehicles  (7  states  each)  plus  two  states  for  user 
equipment  error  sources  [30]. 
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A. 4  DGPS  Error  States 


Table  A. 8  lists  the  DGPS  error  states  developed  by  Negast  [30].  A  total  of  22 
states  are  included  to  model  the  error  characteristics  of  4  space  vehicles  (5  states 
each)  plus  user  equipment  (two  states)  error  sources. 

A. 5  Simulation  Filter  States 

Table  A. 9  lists  the  13  states  used  in  the  PLS  filter  model. 
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Table  A.l  Litton  93-state  INS  Model:  INS  States  1  ^  29 


State 

Number 

nil 

Definition 

1 

SOx 

X-component  of  vector  angle  from  true  to  computer  frame 

2 

■H 

Y-component  of  vector  angle  from  true  to  computer  frame 

3 

80, 

Z-component  of  vector  angle  from  true  to  computer  frame 

4 

X-component  of  vector  angle  from  true  to  platform  frame 

5 

Y-component  of  vector  angle  from  true  to  platform  frame 

6 

Z-component  of  vector  angle  from  true  to  platform  frame 

7 

sv. 

X-component  of  error  in  computed  velocity 

8 

SVy 

Y-component  of  error  in  computed  velocity 

9 

sv. 

Z-component  of  error  in  computed  velocity 

10 

Sh 

Error  in  vehicle  altitude  above  reference  ellipsoid 

11 

ShL 

Error  in  lagged  inertial  altitude 

12 

6S3 

Error  in  vertical  channel  aiding  state 

13 

6S4 

Error  in  vertical  channel  aiding  state 

14 

Oxc 

X-component  of  gyro  correlated  drift  rate 

15 

Oyc 

Y-component  of  gyro  correlated  drift  rate 

16 

bz. 

Z-component  of  gyro  correlated  drift  rate 

17 

^  Xc 

X-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

18 

V,. 

Y-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

19 

V.. 

Z-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

20 

8gx 

X-component  of  gravity  vector  errors 

21 

HOI 

Y-component  of  gravity  vector  errors 

22 

8gz 

Z-component  of  gravity  vector  errors 

23 

Ohs 

Total  baro-altimeter  correlated  error 

24 

bxt 

X-component  of  gyro  trend 

25 

byt 

Y-component  of  gyro  trend 

26 

bzt 

Z-component  of  gyro  trend 

27 

X-component  of  accelerometer  trend 

28 

^yt 

Y-component  of  accelerometer  trend 

29 

Z-component  of  accelerometer  trend 
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Table  A. 2  Litton  93-state  INS  Model:  INS  States  30  —y  47 


Definition 


X-component  of  gyro  drift  rate  repeatability 


Y-component  of  gyro  drift  rate  repeatability 


Z-component  of  gyro  drift  rate  repeatability 


X-component  of  gyro  scale  factor  error 


Y-component  of  gyro  scale  factor  error 


Z-component  of  gyro  scale  factor  error 


X  gyro  misalignment  about  Y-axis 


Y  gyro  misalignment  about  X-axis 


Z  gyro  misalignment  about  X-axis 


X  gyro  misalignment  about  Z-axis 


Y  gyro  misalignment  about  Z-axis 


Z  gyro  misalignment  about  Y-axis 


X  gyro  scale  factor  nonlinearity 


Y  gyro  scale  factor  nonlinearity 


Z  gyro  scale  factor  nonlinearity 


X  gyro  scale  factor  asymmetry  error 


Y  gyro  scale  factor  asymmetry  error 


Z  gyro  scale  factor  asymmetry  error 
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Table  A. 3  Litton  93-state  INS  Model:  INS  States  48  — 69 


State 

Number 


State 

Symbol 


Definition 


48 


V 


X-component  of  accelerometer  bias  repeatability 


49 


V. 


Y-component  of  accelerometer  bias  repeatability 


50 


V, 


Z-component  of  accelerometer  bias  repeatability 


51 


X-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


52 


Y-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


53 


Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


54 


X-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


55 


SqAy 


Y-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


56 


Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


57 


A 


Coefficient  of  error  proportional  to  square 
of  measured  acceleration 


58 


ft 


yy 


Coefficient  of  error  proportional  to  square 


of  measured  acceleration 

59 

fzz 

Coefficient  of  error  proportional  to  square 
of  measured  acceleration 

60 

fxy 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

61 

fxz 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

62 

fyx 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

63 

fyz 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

64 

fzx 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

65 

fzy 

Coefficient  of  error  proportional  to  products  of  acceleration 
along  and  orthogonal  to  accelerometer  sensitive  axis 

66 

X  accelerometer  misalignment  about  Z-axis 

67 

/^2 

Y  accelerometer  misalignment  about  Z-axis 

68 

Z  accelerometer  misalignment  about  Y-axis 

69 

0-3 

Z- accelerometer  misalignment  about  X-axis 

Table  A. 4  Litton  93-state  INS  Model:  INS  States  70  — >  93 


State 

Number 

State 

Symbol 

Definition 

70 

Xq 

X- component  of  accelerometer  bias 
thermal  transient 

71 

V., 

Y-component  of  accelerometer  bias 
thermal  transient 

72 

V., 

Z-component  of  accelerometer  bias 
thermal  transient 

73 

^Xq 

X-component  of  initial  gyro  drift  rate 
bias  thermal  transient 

74 

^yq 

Y-component  of  initial  gyro  drift  rate 
bias  thermal  transient 

75 

bzq 

Z-component  of  initial  gyro  drift  rate 
bias  thermal  transient 

76 

F 

xyz 

X  gyro  compliance  term 

77 

F 

1  xyy 

X  gyro  compliance  term 

78 

F 

J-  xyx 

X  gyro  compliance  term 

79 

F 

±  xzy 

X  gyro  compliance  term 

80 

P 

xzz 

X  gyro  compliance  term 

81 

F 

^  xzx 

X  gyro  compliance  term 

82 

F 

yzx 

Y  gyro  compliance  term 

83 

F 

yzz 

Y  gyro  compliance  term 

84 

F 

■t  yzy 

Y  gyro  compliance  term 

85 

Fyxz 

Y  gyro  compliance  term 

86 

F 

yxx 

Y  gyro  compliance  term 

87 

F 

J-  yxy 

Y  gyro  compliance  term 

88 

F 

J-  zxy 

Z  gyro  compliance  term 

89 

F 

zxx 

Z  gyro  compliance  term 

90 

F 

zxz 

Z  gyro  compliance  term 

91 

F 

^  zyx 

Z  gyro  compliance  term 

92 

F 

zyy 

Z  gyro  compliance  term 

93 

F 

^  zyz 

Z  gyro  compliance  term 
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Table  A. 5  Reduced- Order  INS  System  Model:  INS  States  1  ^  20 


State 

Number 

IIEES3I 

Definition 

LN-93 

State 

PLS 

State 

1 

S9j; 

X-component  of  vector  angle  from  true  to  computer  frame 

1 

1 

2 

SOy 

Y-component  of  vector  angle  from  true  to  computer  frame 

2 

2 

3 

86, 

Z-component  of  vector  angle  from  true  to  computer  frame 

3 

3 

4 

X-component  of  vector  angle  from  true  to  platform  frame 

4 

4 

5 

4'y 

Y-component  of  vector  angle  from  true  to  platform  frame 

5 

5 

6 

<l>z 

Z-component  of  vector  angle  from  true  to  platform  frame 

6 

6 

7 

6V, 

X-component  of  error  in  computed  velocity 

7 

7 

8 

SVy 

Y-component  of  error  in  computed  velocity 

8 

8 

9 

sv. 

Z-component  of  error  in  computed  velocity 

9 

9 

10 

6h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

10 

10 

11 

She 

Total  baro-altimeter  correlated  error 

23 

11 

12 

Shi 

Error  in  lagged  inertial  altitude 

11 

14 

13 

SSs 

Error  in  vertical  channel  aiding  state 

12 

15 

14 

6Si 

Error  in  vertical  channel  aiding  state 

13 

16 

15 

V.. 

X-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

17 

17 

16 

V,. 

Y-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

18 

18 

17 

V.. 

Z-component  of  accelerometer  and 
velocity  quantizer  correlated  noise 

19 

19 

18 

8gx 

X-component  of  gravity  vector  errors 

20 

20 

19 

6gy 

Y-component  of  gravity  vector  errors 

21 

21 

20 

8gz 

Z-component  of  gravity  vector  errors 

22 

22 
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Table  A. 6  Reduced- Order  INS  System  Model:  INS  States  21 


State  Sta 
Number  Sym 


Definition 


X-component  of  gyro  drift  rate  repeatability 


Y-component  of  gyro  drift  rate  repeatability 


Z-component  of  gyro  drift  rate  repeatability 


X-component  of  gyro  scale  factor  error 


Y-component  of  gyro  scale  factor  error 


Z-component  of  gyro  scale  factor  error 


X-component  of  accelerometer  bias  repeatability 


Y-component  of  accelerometer  bias  repeatability 


Z-component  of  accelerometer  bias  repeatability 


X-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


Y-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  error 


X-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


Y-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


Z-component  of  accelerometer  and  velocity 
quantizer  scale  factor  asymmetry 


X  accelerometer  misalignment  about  Z-axis 


Y  accelerometer  misalignment  about  Z-axis 


Z  accelerometer  misalignment  about  Y-axis 


X- accelerometer  misalignment  about  Y-axis 


Y- accelerometer  misalignment  about  X-axis 


Z-accelerometer  misalignment  about  X-axis 
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Table  A. 7  GPS  Error  States 


State 

Sta 

Number 

Sym 

Definition 

PLS 

State 

User  clock  bias 

12 

User  clock  drift 

13 

SV  1  code  loop  error 

44 

SV  1  tropospheric  error 

45 

SV  1  ionospheric  error 

46 

SV  1  clock  error 

47 

SV  1  x-component  of  position  error 

48 

SV  1  y-component  of  position  error 

49 

SV  1  z-component  of  position  error 

50 

SV  2  code  loop  error 

51 

SV  2  tropospheric  error 

52 

SV  2  ionospheric  error 

53 

SV  2  clock  error 

54 

SV  2  x-component  of  position  error 

55 

SV  2  y-component  of  position  error 

56 

SV  2  z-component  of  position  error 

57 

SV  3  code  loop  error 

58 

SV  3  tropospheric  error 

59 

SV  3  ionospheric  error 

60 

SV  3  clock  error 

61 

SV  3  x-component  of  position  error 

62 

SV  3  y-component  of  position  error 

63 

SV  3  z-component  of  position  error 

64 

SV  4  code  loop  error 

65 

SV  4  tropospheric  error 

66 

SV  4  ionospheric  error 

67 

SV  4  clock  error 

68 

SV  4  x-component  of  position  error 

69 

SV  4  y-component  of  position  error 

70 

SV  4  z-component  of  position  error 

71 
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Table  A. 8  DGPS  Error  States 


State 

Number 

IH 

Definition 

PLS 

State 

1 

^Rudkn 

User  clock  bias 

12 

2 

^Duciku 

User  clock  drift 

13 

3 

S  Rf>pQp^ 

SV  1  tropospheric  error 

45 

4 

^Rion\ 

SV  1  ionospheric  error 

46 

5 

SV\ 

SV  1  x-component  of  position  error 

48 

6 

^Vsvi 

SV  1  y-component  of  position  error 

49 

7 

Sz 

SV  1  z-component  of  position  error 

50 

8 

^  R'tTO'P2 

SV  2  tropospheric  error 

52 

9 

SRion2 

SV  2  ionospheric  error 

53 

10 

^^SV2 

SV  2  x-component  of  position  error 

55 

11 

^ysv2 

SV  2  y-component  of  position  error 

56 

12 

Sz 

^^SV2 

SV  2  z-component  of  position  error 

57 

13 

^Rtro'pz 

SV  3  tropospheric  error 

59 

14 

^Rioui 

SV  3  ionospheric  error 

60 

15 

SV  3  x-component  of  position  error 

62 

16 

^VsVi 

SV  3  y-component  of  position  error 

63 

17 

SV  3  z-component  of  position  error 

64 

18 

^Rtropi 

SV  4  tropospheric  error 

66 

19 

^Rion^ 

SV  4  ionospheric  error 

67 

20 

SV  4  x-component  of  position  error 

69 

21 

^Vsv^ 

SV  4  y-component  of  position  error 

70 

22 

SV  4  z-component  of  position  error 

71 
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Table  A. 9  Simulation  Filter  States: 


State 

Number 

nil 

Definition 

1 

86  X 

X-component  of  vector  angle  from  true  to  computer  frame 

2 

86  y 

Y-component  of  vector  angle  from  true  to  computer  frame 

3 

86, 

Z-component  of  vector  angle  from  true  to  computer  frame 

4 

(j>x 

X-component  of  vector  angle  from  true  to  platform  frame 

5 

4'y 

Y-component  of  vector  angle  from  true  to  platform  frame 

6 

<f>Z 

Z-component  of  vector  angle  from  true  to  platform  frame 

7 

sv. 

X-component  of  error  in  computed  velocity 

8 

SVy 

Y-component  of  error  in  computed  velocity 

9 

sv. 

Z-component  of  error  in  computed  velocity 

10 

8h 

Error  in  vehicle  altitude  above  reference  ellipsoid 

11 

Shs 

Total  baro-altimeter  correlated  error 

12 

^Ruclku 

GPS  User  clock  bias 

13 

SDuciku 

GPS  User  clock  drift 
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Appendix  B.  Dynamics  and  Noise  Matrices 

B.l  Model  Dynamics  Matrices  [16,41,4^] 

The  LN-93  error-state  dynamics  matrix  F  as  provided  by  Litton  is  a  93-by-93 
array  that  contains  a  large  number  of  elements  that  are  identically  zero.  Litton 
partitions  the  F  matrix  into  thirty-six  subarrays  [16]  reflecting  the  logical  divisions 
of  error  sources  discussed  in  Chapter  3. 

The  non-zero  elements  of  the  Litton  model  are  included  in  Tables  B.l  through  B.8. 
The  dynamics  submatrices  for  the  4Tstate  INS  model  may  be  constructed  as  needed 
from  the  LN-93  dynamics  submatrices  and  the  full-  to  reduced-order  state  relations 
given  in  Tables  A.l  through  A. 4.  The  reader  should  note  that  Negast’s  [30]  revised 
baro-altimeter  model  states  are  not  included  in  this  set  of  original  F  matrix  elements 
extracted  from  the  Litton  document  [16]. 

A  notational  convention  is  to  label  elements  of  the  C*,  sensor-to-true,  matrix 
as  Cij  where  i  is  the  row  and  j  is  the  column  in  the  transformation  matrix. 

The  dynamics  matrices  for  the  GPS  and  DGPS  truth  models  are  included  in 
the  body  of  Chapter  3. 

B.2  Process  Noise  Matrices 

The  Litton  document  [16]  includes  a  93-by-93  process  noise  matrix  Q  for  the 
LN-93  error  model.  Like  the  F  matrix,  the  Q  matrix  is  partitioned  into  subarrays 
which  correspond  to  the  error-state  subvectors  discussed  in  Chapter  3.  The  vast 
majority  of  the  elements  in  the  process  noise  matrix  are  identically  zero.  The  non¬ 
zero  elements  of  Q  are  shown  in  Tables  B.9  and  B.IO. 

The  process  noise  matrices  for  the  GPS  and  DGPS  models  are  included  in  the 
body  of  Chapter  3,  while  the  Q  noise  strengths  used  to  tune  the  13-state  DGPS  Alter 
model  for  each  of  the  four  navigation  component  cases  is  given  in  Appendix  C. 
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Table  B.l  Elements  of  the  Dynamics  Submatrix  Fn 


Element 

Term 

Element 

Term 

(1.3) 

~Py 

(1.8) 

—Cry 

(2.3) 

Px 

(2,7) 

Crx 

(3.1) 

Py 

(3,2) 

px 

(4.2) 

-fiz 

(4,3) 

fly 

(4.5) 

^iuz 

(4,6) 

^iriy 

(4.8) 

—Cry 

(5.1) 

flz 

(5.3) 

(5,4) 

-(^inz 

(5.6) 

(5,7) 

Crx 

(6.1) 

—  fly 

(6.2) 

fix 

(6.4) 

<^iny 

(6,5) 

(7.1) 

-2K,Dy  -  2V^Ct^ 

(7,2) 

(7.3) 

2Vzny 

(7,5) 

-A, 

(7,6) 

Ay 

(7.7) 

-VzCrx 

(7.8) 

(7,9) 

Py  2,fly 

(8,1) 

2V^ny 

(8,2) 

-2VA  -  2Vzflz 

(8.3) 

2Vzfly 

(8,4) 

Az 

(8,6) 

(8,7) 

-2flz 

(8,8) 

—  VzCry 

(8,9) 

px  +  20a, 

(9.1) 

2K0, 

(9.2) 

2Vyflz 

(9,3) 

-2Vyfly  -  2140^ 

(9,4) 

—  Ay 

(9,5) 

A,; 

(9,7) 

Py  -1-  2fly  A  VxCrX 

(9,8) 

““ Px  20a;  HY 

(9,10) 

2gola 

(9,11) 

-h 

(9,12) 

-1 

(9,13) 

h 

(10,9) 

1 

(10,11) 

-h 

(10,13) 

ki  —  \ 

(11,10) 

1 

(11,11) 

-1 

(12,11) 

h 

(12,13) 

-h 

(13,10) 

ki 

(13,11) 

—k^ 

(13,13) 

ki  —  l 

px,y  =  Components  of  angular  rate,  nav  reference  frame  to  earth-fixed  frame 

^x,y,z  =  Components  of  angular  rate,  earth-fixed  frame  to  inertial  frame 

y  z  —  Components  of  angular  rate,  nav  reference  frame  to  inertial  frame 
14, 2/, z  =  Components  of  vehicle  velocity  vector  in  earth-fixed  coordinates 

■'^x,y,z  =  Components  of  specific  force  in  the  sensor  reference  frame 

^1,2, 3, 4  =  Vertical  channel  gains.  See  LN-93  documentation  [17]  for  equations 

Crx,ry  =  Components  of  earth  spheroid  inverse  radii  of  curvature 
go  =  Equatorial  gravity  magnitude  (32.08744/t/sec^) 

a  =  Equatorial  radius  of  the  earth  (6378388m) 
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Table  B.2  Elements  of  the  Dynamics  Submatrix  F12 


Element 

Term 

Element 

Term 

Element 

Term 

(4,14) 

Cn 

(4,15) 

C12 

(4,16) 

Cl3 

(4,24) 

Cut 

(4,25) 

Cl2t 

(4,26) 

Cist 

(5,14) 

C21 

(5,15) 

C22 

(5,16) 

C23 

(5,24) 

C2lt 

(5,25) 

C22t 

(5,26) 

C2st 

(6,14) 

C3I 

(6,15) 

C32 

(6,16) 

Css 

(6,24) 

Csit 

(6,25) 

C32t 

(6,26) 

Csst 

(7,17) 

Cn 

(7,18) 

C12 

inma 

Cis 

(7,20) 

1 

(7,27) 

Cut 

in^B 

Cl2t 

(7,29) 

Cist 

(8,17) 

C21 

(8,18) 

C22 

(8,19) 

C23 

(8,21) 

1 

(8,27) 

C2lt 

(8,28) 

C22t 

(8,29) 

C23t 

(9,17) 

Csi 

(9,18) 

C32 

(9,19) 

C33 

(9,22) 

1 

(9,23) 

k2 

(9,27) 

Czit 

(9,28) 

Cs2t 

(9,29) 

Czzt 

(10,23) 

h 

(12,23) 

-ks 

(13,23) 

him 

Note:  For  the  above  element  definitions  to  =  0 
C  =  Coordinate  transformation  matrix,  body  frame  to  nav  reference  frame 
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Table  B.3  Elements  of  the  Dynamics 


Element 

Term 

Element 

Term 

(4,30) 

Cn 

(4,31) 

Cl  2 

(4,33) 

(4.34) 

Cui^iny 

(4,36) 

CixCOiriz 

(4.37) 

—  C\2^inz 

(4,39) 

(4,40) 

Cl2<^inx 

(4,42) 

(4,43) 

(4,45) 

(4,46) 

(5,30) 

Car 

(5,31) 

C22 

(5,33) 

(5.34) 

(5,36) 

(5.37) 

—  C22‘^inz 

(5,39) 

(5,40) 

C22^inx 

(5,42) 

(5,43) 

(5,45) 

(5,46) 

(6,30) 

C'si 

(6,31) 

C32 

(6,33) 

(6,34) 

C32<^iny 

(6,36) 

CsiUin^ 

—C32<^inz 

(6,39) 

IH^im 

(6,42) 

(6,43) 

(6,45) 

O.5C31  <^inj,\ 

(6,46) 

0.5C32|a;i„^  | 
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Submatrix  F13 


Element 


(4,32) 


(4,35) 


(4,38) 


(4,41) 


(4,44) 


(4,47) 


(5,32) 


(5,35) 


(5,38) 


(5,41) 


(5,44) 


(5,47) 


C'l3 


^IS^in 


Ci3U>in. 


-ClSUi 


(6,38) 


(6,41) 


(6,44) 


C33 


CssUJiny 


C33U^L 


Table  B.4  Elements  of  the  Dynamics  Submatrix  F14 


Element 


Term 


Element 


Term 


Element 


Term 


(7,48) 


C 


11 


(7,49) 


C 


12 


(7,50) 


C 


13 


(7,51) 


CnAS 


(7,52) 


CioA; 


B 


(7,53) 


(7,54) 


Cu\A^\ 


(7,55) 


Ci2\Ay\ 


(7,56) 


C'lalAf  I 


(7,57) 


CuA^A^ 


(7,58) 


CrsA^Af 


(7,59) 


(7,60) 


X 


(7,61) 


(7,62) 


(7,63) 


CuA^A 


(7,64) 


-C-ioAf 


(7,65) 


a 


(7,66) 


(7,67) 


(7,68) 


(7,69) 


^7i341 


(8,48) 


C 


21 


(8,49) 


0 


22 


C22A^ 


(8,50) 


<^23 


(8,51) 


(8,52) 


(8,53) 


(8,54) 


f72i!Afl 


'W 

X 


(8,55) 


C^22|Af| 


(8,56) 


C23U 


iB'i 


(8,57) 


<^2iA 


(8,58) 


C22A 


(8,59) 


C^zAf 
C22A^A^ 


(8,60) 


C2l4lfAf 


C'22A^Af 


(8,61) 


C2iA^A? 


(8,62) 


^23^: 


(8,63) 


(8,64) 


(8,65) 


(8,66) 


¥ 


(8,67) 


-C22A^ 


(8,68) 


(8,69) 


C23A 


(9,48) 


Csi 

C31A" 


(9,49) 


a 


32 


(9,50) 


C; 


33 


(9,51) 


(9,52) 


a324f_ 


(9,53) 


C33A?' 


(9,54) 


^3l|Af| 


W 

X 


(9,55) 


C32|Af  I 


(9,56) 


C'33|Af'| 


(9,57) 


CsiA: 


(9,58) 


C32A?' 


(9,59) 


C33A 


Cs2A^A^ 


(9,60) 


(9,61) 


C: 


31 


CssA^A^ 
—  C32A^ 


(9,62) 


C33A^A^ 


(9,63) 


a 


32 


(9,64) 


(9,65) 


^'33^4?^ 


(9,66) 


¥ 


(9,67) 


(9,68) 


(9,69) 


C'asA 


^^,y,z  —  Components  of  acceleration  in  the  body  frame 

Af '  =  Specific  force  component  (includes  gravity) 
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Table  B.5  Elements  of  the  Dynamics  Submatrix  F15 


Element 

Term 

Element 

Term 

Element 

Term 

(4,73) 

Cn 

(4,74) 

C\2 

(4,75) 

Cl3 

(6,73) 

C21 

(5,74) 

C22 

(5,75) 

C23 

(6,73) 

C-si 

(6,74) 

Cz2 

(6,75) 

^33 

(7,70) 

Cn 

(7,71) 

C\2 

(7,72) 

(713 

(8,70) 

C21 

(8,71) 

C22 

(8,72) 

C23 

(9,70) 

C31 

(9.71) 

C32 

(9,72) 

C33 

Table  B.6  Elements  of  the  Dynamics  Submatrix  Fie 


Element 

Term 

Element 

Term 

Element 

Term 

(4,76) 

(4,77) 

(4,78) 

ISBEBfflB 

(4,79) 

(4,80) 

(4,81) 

^^BM^BB 

(4,82) 

(4,83) 

(4,84) 

^^BM^BB 

(4,85) 

(4,86) 

(4,87) 

(4,88) 

(4,89) 

(4,90) 

IgffijilaflBM 

(4,91) 

(4,92) 

(4,93) 

(5,76) 

(5,77) 

(5,78) 

(5,79) 

(5,80) 

(5,81) 

(5,82) 

(5,83) 

(5,84) 

KMBtaflBBI 

(5,85) 

(5,86) 

(5,87) 

(5,88) 

(5,89) 

(5,90) 

(5,91) 

(5,92) 

(5,93) 

BEBCTafl^M 

(6,76) 

(6,78) 

IgRP^yB^M 

(6,79) 

iimwi 

(6,81) 

(6,82) 

(6,83) 

(6,84) 

(6,85) 

C 32  A^Uin^ 

(6,86) 

(6,87) 

(6,88) 

(6,89) 

(6,90) 

IggHiUW 

(6,91) 

(6,92) 

(6,93) 
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Table  B.7  Elements  of  the  Dynamics  Submatrix  F22 


Element 

Term 

Element 

Term 

Element 

Term 

(14,14) 

(15,15) 

(16,16) 

-l^bxr. 

(17,17) 

(18,18) 

(19,19) 

-I^^zr. 

(20,20) 

~^Sgx 

(21,21) 

~I^Sgy 

(22,22) 

-I^Sgx 

(23,23) 

-^Shc 

/^V 

xc,yc,zc 

I^Shc 


=  Gyro  inverse  correlation  time  constants  (5  min) 

=  Accelerometer  inverse  correlation  time  constants  (5  min) 

=  Gravity  vector  error  inverse  correlation  time  constants  (V /20NM) 
=  Barometer  inverse  correlation  time  (10  min) 


Table  B.8  Elements  of  the  Dynamics  Submatrix  F55 


Element 

Term 

Element 

Term 

Element 

Term 

(70,70) 

-P^xa 

(71,71) 

(72,72) 

(73,73) 

-^bxa 

(74,74) 

-^b,„ 

(75,75) 

By  =  Accelerometer  bias  thermal  transient  inverse  time  constants  (1  min) 

n  V  xq,yq,zq  \  / 

Bb:cq.yq.zq  =  Gyro  drift  rate  bias  thermal  transient  inverse  time  constants  (1  min) 
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Table  B.9  Non-zero  Elements  of  Process  Noise  Submatrix  Qn 


Element 

Term 

Element 

Term 

(4,4) 

(5,5) 

(6,6) 

<- 

(7,7) 

< 

(8,8) 

(9,9) 

< 

Table  B.IO  Non-zero  Elements  of  Process  Noise  Submatrix  Q22 


Element 

Term 

Element 

Term 

(14,14) 

2^6.. 

(15,15) 

(16,16) 

(17,17) 

2/5v.,<Tv., 

(18,18) 

(18,18) 

(20,20) 

(21,21) 

(22,22) 

‘^I^SgySa, 

(23,23) 

^bx 


VAx 


a 


°Sx,y,z 

2 

She 


=  PSD  value  of  gyro  drift  rate  white  noise  (6.25e-10^^) 

=  PSD  value  of  accelerometer  white  noise  (1.037e-7^-) 

=  Variances  of  gyro  drift  correlated  noise  (3.086e-13^^) 

=  Variances  of  accelerometer  correlated  noise  (4.147e-9-^) 

=  Variances  of  gravity  vector  error  component  correlated  noise  (1.93e-6de£f^ 
=  Variance  of  barometer  correlated  noise  (10000/i^) 
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Appendix  C.  Filter  Tuning  Values 


This  appendix  contains  a  listing  of  the  dynamics  driving  noise  strength  val¬ 
ues  used  to  tune  the  Kalman  filters  for  each  navigation  component  case  shown  in 
Table  C.l.  “Tuning”,  or  purposeful  addition  of  dynamics  driving  noise,  is  done  to 
account  for  the  modelling  inadequacies  introduced  by  using  a  small  number  of  states 
(13  in  this  study)  to  represent  a  high-order  truth  model  (95  or  96  states  in  this 
study),  which  in  turn  uses  (usually  much  smaller)  values  of  dynamics  driving  noise 
strengths  to  account  for  its  reduction  from  the  infinite-state  real-world  model.  The 
95-state  truth  model  developed  by  Litton  for  the  LN-93  INS  [17]  is  accompanied  by 
the  dynamics  driving  noise  strength  values  shown  in  the  variable  definitions  below 
Table  B.IO. 

For  the  13  states  in  the  filter  model  used  in  this  application  (see  Table  A. 9), 
tuning  was  accomplished  by  experimentally  adjusting  the  Q  values  shown  in  Ta¬ 
ble  C.2  for  each  of  the  four  navigation  component  cases.  The  13  Q  rows  of  Table  C.2 
correspond  directly  to  the  13  states  of  the  filter  model  as  defined  in  Table  A. 9. 


Table  C.l  Navigation  Component  Cases 


Case 

INS  Type 

GPS  Type 

Altimeter  Aiding 

1 

0.4nm/hr 

DGPS  and 
Pseudolite 

Baro  and 
Radar  Alt 

2 

0.4nm/hr 

DGPS 

Baro  and 
Radar  Alt 

3 

0.4nm/hr 

DGPS 

Baro 

4 

4.0nm/hr 

DGPS 

Baro 
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Table  C.2  Tuning  Values  for  Filter  States 


State  No. 

Nav  Case  1 
Values 

Nav  Case  2 
Values 

Nav  Case  3 
Values 

Nav  Case  4 
Values 

1 

1.2e-16 

4.0e-15 

4.0e-16 

1.2e-23 

2 

1.5e-16 

3.0e-15 

3.0e-16 

1.5e-23 

3 

0.0 

0.0 

0.0 

0.0 

4 

3.81e-12 

1.43e-12 

1.43e-12 

5.71e-9 

5 

5.33e-12 

3.81e-12 

3.81e-12 

5.71e-9 

6 

1.62e-ll 

2.38e-12 

2.38e-12 

3.43e-10 

7 

1.03e-8 

1.29e-5 

1.29e-5 

5.15e-5 

8 

1.29e-5 

1.29e-5 

1.03e-8 

9 

2.78e-3 

2.78e-3 

2.98e-3 

10 

16.0 

15.0 

24.0 

11 

6.67e3 

6.67e3 

3.33e4 

12 

0.2 

7.5 

7.5 

0.2 

13 

5.0e-15 

S.Oe-lO 

5.0e-10 

5.0e-16 
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